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ABSTRACT 


CHASING  SHADOWS  IN  THE  OUTER  SOLAR  SYSTEM 


Federica  Bianco 


Charles  R.  Alcock 


The  characteristics  of  the  populations  of  objects  that  inhabit  the  outer 
solar  system  carry  the  hngerprint  of  the  processes  that  governed  the  formation  and 
evolution  of  the  solar  system.  Occultation  surveys  push  the  limit  of  observation 
into  the  very  small  and  distant  outer  solar  system  objects,  allowing  us  to  set  con¬ 
straints  on  the  structure  of  the  Kuiper  belt,  Scattered  disk  and  Sedna  populations. 
I  collected,  reduced,  and  analyzed  vast  datasets  looking  for  occultations  of  stars  by 
outer  solar  system  objects,  both  working  with  the  Taiwanese  American  Occultation 
Survey  (TAOS)  collaboration  and  leading  the  MMT/Megacam  occultation  effort. 
Having  found  no  such  events  in  my  data,  I  was  able  to  place  upper  limits  on  the 
Kuiper  belt,  scattered  disk  and  Sedna  population.  These  limits  and  their  derivation 
are  described  here. 


V 


Contents 


List  of  Tables  viii 

List  of  Figures  ix 

1  Introduction  1 

1.1  Historical  note .  1 

1.2  Formation .  2 

1.3  Observational  techniques  to  explore  the  outer  solar  system .  6 

1.3.1  Diffraction  dominated  occultations  of  bright  stars .  8 

1.4  The  story  told  by  the  small  KBOs .  15 

2  The  sub-km  end  of  the  Kuiper  Belt  size  distribution  20 

2.1  Introduction . 21 

2.2  Fast  Photometry  with  a  Large  telescope:  The  Continuous  -  Readout 

Mode . 23 

2.3  Data . 26 

2.3.1  Data  extraction  and  reduction . 28 

2.4  Residual  noise  in  the  time-series . 35 

2.5  Search  for  events  and  efficiency  . 40 

2.5.1  Detection  algorithm . 40 

2.5.2  Efficiency . 42 

2.5.3  Rejection  of  false  positives . 45 

2.6  Upper  limit  to  the  size  distribution  of  KBOs  and  scientihc  interpretation  46 

2.6.1  Comparison  with  the  results  from  the  TAOS  survey . 48 

2.6.2  The  Kuiper  belt  as  reservoir  of  Jupiter  Family  Comets  ....  49 

3  The  TAOS  survey  and  the  3.75-year  dataset  51 

3.1  Introduction . 51 

3.2  3.75  years  of  TAOS  data . 52 

3.3  Analysis . 54 

3.3.1  Photometry . 57 

3.3.2  Detection  and  false  positive  rejection . 58 


VI 


3.3.3  Occultation  event  simulator . 63 

3.3.4  Determination  of  the  stellar  angular  size . 69 

3.3.5  Implantation  and  efficiency  test  parameters  . 71 

3.3.6  Analysis  of  the  efficiency  parameters . 73 

4  Constraints  on  models  of  the  Solar  System  formation  and  evolntion 

from  the  TAOS  data  77 

4.1  Introduction . 78 

4.2  Effective  coverage  and  upper  limits . 78 

4.3  The  Jupiter  Family  Comets  progenitor  population . 79 

4.4  Outer  Solar  System  collisional  models . 82 

4.4.1  Pan  &  Sari  (2005)  .  84 

4.4.2  Kenyon  &  Bromley  (2004) .  86 

4.4.3  Benavidez  &  Campo  Bagatin  (2009)  .  88 

4.4.4  Generic  3-regime  model:  constraints  on  the  intermediate  re¬ 
gion  of  the  size  spectrum . 89 

5  Exploring  the  Solar  System  beyond  the  Kuiper  belt  94 

5.1  Introduction . 94 

5.2  Search  algorithm . 98 

5.3  Event  rate  calculation  for  Sedna-like  objects . 100 

5.4  Renewed  limits  to  the  Sedna  population  . 103 

5.5  Future  work . 105 

6  Conclusions  106 

6.1  The  MMT/Megacam  Survey  for  small  KBOs . 106 

6.2  The  TAOS  KBO  survey  . 108 

6.3  The  TAOS  Sedna  and  Scattered  Disk  Survey . 110 


vii 


List  of  Tables 


2.1  MMT/Megacam  survey  observed  fields . 27 

2.2  MMT/Megacam  survey  data  set  parameters . 28 

3.1  TAOS  dataset  parameters  (3-telescope  data) . 55 

3.2  Distribution  of  synthetic  events . 72 

3.3  Parameter  of  implanted  events  in  Figure  3.7 . 74 


viii 


List  of  Figures 


1.1  Outer  Solar  System  taxonomy .  6 

1.2  TNO  discoveries  by  year  and  R  magnitude .  8 

1.3  Luminosity  distribution  of  KBOs  from  Fuentes  &  Holman  (2008)  ...  9 

1.4  Occultation  geometry .  11 

1.5  Diffraction  by  a  TNO  ocultation .  13 

1.6  Diameter  -  semimajor- axis  phase  space  detectability  . 15 

2.1  Simulated  diffraction  pattern  for  an  MMT/Megacam  hypothetical 

target  . 21 

2.2  Megacam  focal  plane . 22 

2.3  Images  generated  by  conventional  camera  use  and  by  continuous- 

readout  . 24 

2.4  Frequency  analysis  of  MMT/Megacam  continuous-readout  data  ...  26 

2.5  Image  motion  and  PSF  width  fluctuation  for  our  continuous-readout 

data . 31 

2.6  MMT/Megacam  continuous-readout  lightcurves  before  and  after  de¬ 
trending  . 34 

2.7  Implantation  of  occultation  signatures  in  our  data . 35 

2.8  Non-parallel,  simultaneous  CCD  motion  features . 39 

2.9  Detection  efficiency  and  number  of  false  positives . 40 

2.10  Detection  efficiency  and  effective  solid  angle  of  our  survey . 45 

2.11  Duration  and  flux  drop  caused  by  occultations  in  the  size  and  distance 

regime  we  are  considering  . 47 

2.12  Upper  limits  to  the  surface  density  of  KBOs  from  the  MMT /Megacam 

occultation  survey . 48 

3.1  Distribution  of  magnitudes  and  SNR  for  the  TAOS  target  stars  ...  55 

3.2  SNR  versus  TAOS  instrumental  magnitude . 56 

3.3  Distribution  of  angles  from  opposition  for  the  TAOS  targets . 56 

3.4  Diffraction  table  and  synthetic  diffraction  features . 65 

3.5  Steps  of  the  generation  of  a  simulated  occultation  event . 68 

3.6  Angular  size  distribution  for  a  typical  TAOS  held . 71 


IX 


3.7  Implanted  occultations  recovered  by  our  pipeline . 73 

3.8  Our  recovery  efficiency  for  3  km  KBOs  as  a  function  of  observational 

parameters . 76 

4.1  Effective  solid  angle  of  our  survey . 79 

4.2  Model  independent  upper  limits  to  the  surface  density  of  KBOs  from 

the  TAOS  survey . 80 

4.3  PS05  model  and  expected  event  rate  for  TAOS . 85 

4.4  Modeling  of  the  KB04  results:  the  D  =  0.5  —  30  km  region . 86 

4.5  KB04  model  and  expected  event  rate  for  TAOS  with  excess  centered 

at  5.5  km  . 91 

4.6  KB04  model  and  expected  event  rate  for  TAOS  with  excess  centered 

at  1.6  km  . 92 

4.7  BCB09  model  and  expected  event  rate  for  TAOS  . 93 

4.8  Generic  three-slopes  parametric  modeling  of  the  KBO  size  distribu¬ 
tion  and  constraints  from  the  TAOS  survey . 93 

5.1  Occultation  signatures  of  a  77  =  10  km  and  a  77  =  5  km  objects  at 

400  AU  and  their  detectability . 100 

5.2  Minimum  coverage  of  our  survey . 102 

5.3  Maximum  number  density  of  77  >  1  km  objects  allowed  by  our 

survey  between  100  and  1000  AU . 104 


5.4  Limits  to  the  density  and  absolute  number  of  Sednas  from  our  survey  105 


X 


Chapter  1 


Introduction 

1.1  Historical  note 

It  could  be  argued  that  interest  in  the  solar  system  initiated  scientihc 
thinking.  All  ancient  cultures  that  left  a  written  or  graphical  record  left  some  rep¬ 
resentation  of  the  Sun  and  the  Moon,  of  the  planets  and  the  hxed  stars.  Dynamical 
theories  of  what  was  then  thought  of  as  the  entire  Universe  are  among  the  hrst 
records  of  most  cultures.  Early  Americans,  Greeks,  Egyptians,  Chinese,  Babylo¬ 
nians,  Indian,  Celtic,  and  Islamic  cultures,  all  generated  cosmological  theories  to 
describe  and  explain  the  rising  and  setting  of  the  sun,  the  phases  of  the  moon,  the 
motions  of  the  planets  across  the  sky,  the  immobile  stars  and  the  changing  of  sea¬ 
sons.  The  motions  of  solar  system  objects,  which  we  now  call  ephemerides,  were 
studied  well  enough  by  the  ancient  Egyptians  that  they  were  able  to  build  enormous 
tomb  structures  that  allow  the  Sun  to  shine  in  precisely  once  a  year  and  only  once 
a  year.  Mechanical  tools  (e.g.:  the  Antikythera)  were  designed  to  simplify  the  cal¬ 
culation  of  an  ephemeris,  and  are  considered  today  to  be  precursors  to  the  modern 
calculators  and  computers  (Write,  2007,  and  references  therein). 
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Ancient  Greek  philosophers  payed  great  attention  to  astronomy,  and  pro¬ 
duced  a  variety  of  theories  to  describe  the  motion  of  objects  in  the  sky,  including 
some  early  heliocentric  theories  (Aristarchus  of  Samos,  Batten  1981),  untill  the  Aris¬ 
totelian  idea  of  concentric  spheres  prevailed  and  dominated  in  some  form  or  another, 
until  Copernicus  and  Galileo^.  Galileo,  Kepler,  Newton,  and  later  Kant,  started  in¬ 
vestigating  the  forces  which  keep  planets  moving  in  their  orbits,  what  holds  them 
up  in  place  and  puts  them  in  motion. 

Remarkably,  after  over  2000  years  of  solar  system  science,  we  have  probed 
only  as  little  as  10“®  of  the  space  occupied  by  the  solar  system:  the  inner  portion, 
which  contains  over  99%  of  its  mass.  The  outskirts  of  the  solar  system  are  however 
still  unexplored.  Outside  of  the  orbit  of  Neptune  the  solar  system  is  populated 
by  small  (thousands  of  kilometers  down  to  dust  grain  size)  icy  bodies,  arranged  in 
different  structures:  the  Kuiper  belt,  the  Oort  cloud,  and  the  scattered  disk.  Our 
direct  observational  knowledge  is  limited  to  relatively  large  objects  populating  the 
Kuiper  belt  and  the  scattered  disk. 

1.2  Formation 

In  the  current  solar  system  formation  and  evolution  scenario  planets  formed 
by  gravitational  instabilities  in  the  protoplanetary  nebula:  a  disk  of  gas  and  dust 
surrounding  the  Sun. 

Throughout  the  remainder  of  this  section  I  will  refer  to  the  solar  system 
evolution  model  that  came  to  be  known  as  the  Nice  model  (Gomes  et  ah,  2005a; 
Morbidelli  et  ah,  2005;  Tsiganis  et  ah,  2005;  Levison  et  ah,  2008,  and  references 
therein),  one  of  the  most  successful  scenarios  to  explain  the  current  conhguration  of 

^Aristotle  (Greek  philosopher)  -  Britannica  Online  Encyclopedia,  Britannica.com. 
http://www.britannica.com/EBchecked/topic/34560/Aristotle.  Retrieved  March  2,  2010. 
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solar  system  objects.  In  this  model  the  four  giant  planets,  Jupiter,  Saturn,  Uranus, 
and  Neptune,  originally  formed  from  the  circumsolar  disk  of  dust  and  gas  on  nearly 
circular  orbits  between  5  and  17  AU,  a  much  more  compact  conhguration  than  that 
which  we  see  now.  Gravitational  instabilities  in  the  disk  caused  accretion  processes 
to  form  large  planets  out  of  the  dust  and  small  planetesimals  inhabiting  the  disk. 
Outside  of  the  planetary  orbits  there  remained  a  disk  of  icy  planetesimals,  extending 
out  to  ~  35  AU. 

In  the  early  chaotic  solar  system  the  planetesimals  which  were  residing 
within  the  planetary  orbits  suffered  occasional  gravitational  encounters  with  the 
outermost  giant  planets,  which  changed  the  planetesimal  orbits,  initially  scattering 
them  inward.  In  turn,  the  outermost  giant  planets  began  migrating  outward  to 
preserve  angular  momentum.  As  they  were  pushed  into  the  inner  regions  of  the  solar 
system,  the  planetesimals  began  to  interact  with  Jupiter,  which  on  account  of  its 
much  larger  mass  dramatically  affected  their  orbital  parameters  at  each  encounter. 
The  icy  bodies  that  encounter  Jupiter  were  sent  onto  highly  elliptical  orbits  or 
ejected  altogether  from  the  solar  system,  and  to  compensate  for  the  loss  of  angular 
momentum  Jupiter  migrated  inward. 

After  several  hundred  million  years  Jupiter  and  Saturn  crossed  their  1:2 
mean  motion  resonance  and  that  suddenly  increased  their  orbital  eccentricities, 
destabilizing  the  entire  solar  system.  Uranus’s  and  Neptune’s  eccentricities  were 
increased  as  a  result  of  the  interaction  with  Saturn,  and  as  they  migrated  outward 
they  pushed  planetesimals  outward  as  well,  wiping  the  region  which  they  transited 
clean  of  icy  bodies.  Through  these  interactions,  the  orbits  of  Uranus  and  Neptune 
were  again  circularized.  Meanwhile  the  vast  majority  of  planetesimals  were  shuffled 
and  over  90%  of  the  original  planetesimals  mass  was  removed. 

Of  the  planetesimals  that  have  been  pushed  outward  by  the  migration  of 
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Neptune,  a  large  portion  remained  in  low  inclination  orbits,  and  some  of  them  got 
captured  in  mean  motion  resonances  with  Neptune.  These  formed  the  Kuiper  belt, 
hrst  predicted  by  Kenneth  Edgeworth  in  1943,  and  later  by  Gerard  Kuiper  in  1951. 
The  Kuiper  belt  is  a  flat  structure,  extending  from  the  orbit  of  Neptune  (~  30  AU) 
out  to  about  50  AU,  where  observations  reveal  an  outer  edge  (Jewitt  et  ah,  1998). 

A  fraction  of  the  Kuiper  belt  planetesimals  were  excited,  possibly  through 
resonance  and  weak  chaos  associated  to  secular  Kozai  mechanisms  (Volk  &  Malhotra, 
2009),  and  their  inclination  distribution  increased,  forming  the  scattered  disk  and 
extended  scattered  disk.  The  scattered  disk  is  a  family  of  TNOs  populated  by 
objects  with  perihelia  g  >  30  AU  and  with  inclinations  as  high  as  i  ~  40°.  The 
extended  scattered  disk  is  a  population  of  object  whose  perihelia  is  large  enough  to 
be  decoupled  from  Neptune  (see  Figure  1.1). 

Many  planetesimals  bounced  in  and  out  of  the  inner  solar  system  by  re¬ 
peated  encounters  with  Jupiter.  These  interactions  increased  the  eccentricities  e 
of  their  orbits.  The  planetesimals  that  were  not  ejected  via  this  mechanism  were 
hand  off  out  to  Neptune.  Interactions  with  Neptune  continued  increasing  the  or¬ 
bital  sami-major  axis  while  preserving  the  perihelion  distance  of  these  planetesimals. 
Once  their  orbits  reached  ~  3,  000  AU  the  influence  of  the  mean  gravitational  held 
of  the  Galaxy  became  important,  and  their  perihelia  lifted  from  the  region  of  inhu- 
ence  of  Neptune.  The  inclination  distribution  of  planetesimals  in  orbit  farther  than 
3,  000  AU  increased  dramatically:  these  planetesimals  formed  the  Oort  cloud,  an 
isotropic  structure  of  icy  bodies  which  might  extend  farther  out  than  10^  AU  (Dun¬ 
can  et  ah,  1987). 

The  existence  of  a  cloud  of  icy  bodies  in  the  out  skirts  of  the  solar  system 
was  originally  proposed  by  J.  Oort  in  1950  (Oort,  1950),  in  order  to  explain  the  long 
period  the  comets  we  observe.  Long  period  comets  are  indeed  Oort  cloud  objects  at 
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their  perihelion  passage  and  they  provide  the  only  observational  evidence  we  have 
of  the  Oort  cloud. 


Figure  1.1  shows  a  schematic  partition  of  the  perihelion  and  semi- major 
axis  phase  space  and  the  TNO  families  that  occupy  it.  The  classical  Kuiper  belt 
(KB)  lies  on  the  bottom  left  corner  of  the  plot.  An  empty  square  shows  the  position 
of  the  median  perihelion  and  semi-major  axis  of  the  classical  and  resonant  KBOs. 
At  g  >30  AU  and  a  >  50  AU  lays  the  scattered  disk  (SD).  The  extended 

scattered  disk  comprises  objects  with  perihelia  larger  than  g  ~  38  AU,  which  are 
decoupled  from  Neptune.  The  inner  Oort  cloud  start  outward  of  a  =  3,  000  AU 
and  the  outer  Oort  cloud  at  a  >  20,  OOOAU.  The  reader  is  cautioned  that  this 
schematic  separation  in  g  —  a  space  is  only  for  reference.  Other  families  are  here 
ignored,  such  as  the  resonant  populations,  Plutinos,  and  Centaurs.  Furthermore, 
the  TNO  taxonomy  is  not  unique,  and  TNO  nomenclature  is  often  based  on  other 
variable,  g  —  e,  g  —  i,  color,  etc.  (Barucci  et  ah,  2005;  Gladman  et  ah,  2008,  and 
references  therein) 

The  Nice  model  model  is  able  to  explain  many  features  of  our  solar  system; 
it  reproduces  the  current  planetary  orbits,  the  existence  of  objects  in  the  Kuiper 
belt  and  in  resonance  with  Neptune,  the  asteroid  belt,  as  a  relic  of  the  planetary 
nebula  not  affected  by  the  migration  of  Uranus  and  Neptune.  It  is  however  unable 
to  explain  the  existence  of  an  outer  edge  in  the  Kuiper  belt,  it  predicts  more  mass  to 
be  left  in  the  solar  system  than  we  know  of,  and  it  cannot  explain  the  the  existence 
of  Sedna. 

Sedna  was  discovered  in  2004  (Brown  et  ah,  2004).  None  of  the  formation 
mechanisms  in  the  literature  at  the  time  of  discovery  were  able  to  place  orbits  in  this 
region  of  the  perihelion  and  semi-major  axis  phase  space:  this  region  is  inaccessible 
to  objects  coming  from  the  Kuiper  belt  as  well  as  for  objects  being  kicked  out  from 
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Figure  1.1:  Schematic  partition  of  the  perihelion  and  semi-major  axis  phase  space.  Five 
distinct  dynamical  families  are  identified:  the  Kuiper  belt  (KB)  in  the  bottom  left  corner, 
the  scattered  disk  (SD),  the  extended  scattered  disk  (ESD),  the  inner  and  outer  Oort 
cloud.  The  largest  observed  TNOs  are  shown  (black  circles),  as  well  as  the  median  values 
of  q  and  a  of  the  Kuiper  belt  (open  square). 

the  inner  solar  system.  To  date  no  surveys  have  detected  any  other  object  in  orbit 
similar  to  Sedna,  and  its  existence  is  still  unexplained  in  the  current  formation  and 
evolution  scenarios  without  invoking  the  presence  of  a  perturbation  from  an  external 
body  (see  Section  5). 


1.3  Observational  techniques  to  explore  the  outer 
solar  system 

Populated  by  small  and  cold  bodies  the  outer  solar  system  is  among  the 
most  challenging  observational  targets  in  astrophysics  today.  Direct  detection  of 
Trans  Neptunian  Objects  (TNOs)  is  a  difficult  task.  These  objects  typically  range 
from  ~  1000  km  in  diameter,  to  sizes  as  small  as  dust  grains.  They  do  not  shine,  but 
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are  seen  in  reflected  Sun  light,  thus  getting  fainter  as  ~  A^,  where  A  is  the  distance 
to  the  Earth.  Furthermore,  they  move  across  the  sky  at  a  rate  of  ~  3  arcsec/hour, 
making  it  impossible  to  increase  their  detectability  just  by  increasing  the  integration 
time,  and  rendering  any  technique  used  to  increase  the  depth  of  a  survey,  such  as 
stacking,  much  harder  to  perform. 

The  very  Erst  Kuiper  belt  object  to  be  observed  was  Pluto  (134340  Pluto), 
discovered  in  1930  by  Clyde  Tombaugh.  Pluto  has  a  mean  magnitude  R  ~  14, 
and  a  diameter  of  about  2,390  km.  For  its  size,  Pluto  is  exceptionally  bright  due 
to  a  high  albedo  of  about  50%.  Classified  as  a  planet  until  2006,  Pluto  is  today 
the  second  largest  known  KBO  after  Eris  (136199  Eris),  which  has  a  diameter  of 
about  2,500  km  and  an  apparent  magnitude  of  i?  ~  18.4.  Eris  was  discovered  in 
2005  (Brown  et  ah,  2005). 

Over  30  years  passed  between  the  hrst  detection  of  Pluto  and  the  discovery 
of  the  next  KBO.  Jewitt  and  Luu  announced  the  “Discovery  of  the  candidate  Kuiper 
belt  object”  in  a  Nature  paper  on  August  30th  1992,  the  object  known  as  15760 
1992  QBl.  Six  months  later  they  reported  a  second  object  in  the  Kuiper  belt 
region,  181708  1993  FW.  Today  (March  2,  2010)  1099  observed  TNOs  are  cataloged 
by  the  Minor  Planet  Center^,  which  keeps  track  of  all  TNOs  and  minor  planet 
observations  and  creates  ephemerides  to  predict  their  position  in  time.  The  rate  of 
TNO  discoveries,  however,  peaked  in  2001  and  it  has  been  decreasing  ever  since  (see 
Figure  1.2,  left). 

To  date  all  observed  TNOs  are  brighter  than  R  >  30.  Figure  1.2,  right, 
shows  the  /^-magnitude  distribution  of  known  KBOs.  The  majority  of  the  observed 
objects  are  in  the  magnitude  range  23  <  R  <  26;  the  median  magnitude  of  observed 
TNOs  is  {R)  ~  23.5.  Figure  1.3,  from  Fuentes  &  Holman  (2008),  shows  the 

^http://www.cfa. harvard,  edu/iau/ mpc.html 
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Figure  1.2:  Left:  number  of  discovered  TNOs  per  year.  Right:  number  of  observed 
TNOs  as  a  function  of  R  magnitude.  The  vast  majority  of  observed  TNOs  have  apparent 
magnitude  in  the  range  22  <  R  <  26.  Brighter  KBOs  are  rare  and  fainter  KBOs,  although 
more  numerous,  are  hard  to  detect 

cumulative  size  distribution  of  KBOs  as  a  function  of  R  magnitude.  Most  observed 
KBO  are  individually  plotted.  The  density  of  observed  KBOs  with  23  <  i?  <  26 
allows  a  hrm  determination  of  the  size  distribution  in  this  magnitude  range  (see 
Section  1.4).  Only  three  objects  have  been  observed  that  are  fainter  than  magnitude 
R  =  26.5:  R  =  26.7,  28.0  and  28.2,  corresponding  to  diameters  of  respectively  44, 
28  and  25  km  assuming,  as  customary,  an  albedo  of  4%  (Bernstein  et  ah,  2004). 

Future  all  sky  surveys  such  as  Pan-STARRS^  and  LSST"^  will  discover 
many  more  R  <  26  objects  (Jewitt,  2003),  and  data  mining  projects  are  in  progress 
to  detect  faint  targets  in  archival  HST  data  (Fuentes  et  ah,  2009a).  The  observa¬ 
tional  barrier  at  i?  ~  30  is  however  hard  to  overcome. 


1.3.1  Diffraction  dominated  occultations  of  bright  stars 

Bailey  (1976)  proposed  that  small  TNOs  could  be  seen  indirectly  at  their 
passage  across  the  line  of  sight  to  a  star.  This  event,  today  known  as  an  occultation, 
would  produce  a  variation  in  the  flux  of  the  observed  star  which  in  principle  can 


^  http:/ /pan-starrs.ifa.hawaii.  edu/public/science-goals/solar-system.  html 
http://lsst.org/lsst/science/scientist-Outer_solar_system 
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Figure  1.3:  From  Fuentes  &:  Holman  (2008),  this  plot  shows  the  size  distribution  of  large 
KBOs  as  a  function  of  magnitude  (converted  in  diameter  on  the  top  axis).  The  data-points 
include  most  KBOs  observed  directly  to  date. 
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be  observed,  much  like  in  planetary  transit  surveys.  The  dramatic  difference  from 
planetary  transit  surveys  is  however  in  the  geometry  of  the  system;  in  this  technique 
a  distant  star,  often  an  unresolved  point  source,  is  occulted  by  a  near-by  object 
with  a  finite  angular  size.  Noticeably,  for  the  events  of  interest,  the  geometry  of 
the  system  is  such  that  occultations  by  outer  solar  system  objects  too  small  to  be 
observed  directly  are  typically  diffraction  dominated  events. 

If  we  consider  a  roughly  kilometer  size  object  in  the  Kuiper  belt  (30-50  AU) 
we  have  the  case  of  a  light  wave  obstructed  by  an  object  at  finite  distance,  where 
the  diffractor  size  is  large  compared  to  the  wavelength.  In  this  regime  diffraction  is 
properly  described  in  terms  of  the  Huygens-Fresnel  principle.  The  discussion  that 
follows  is  based  on  Born  &  Wolf  (1980)  and  on  Roques  et  al.  (1987). 

Modeling  the  occulting  object  as  a  flat  opaque  screen  S,  and  neglecting 
the  scattering  of  light  at  the  edges  of  the  occulter,  the  diffraction  amplitude  as  of  a 
monochromatic  plane  wave  at  wavelength  A  can  be  derived  at  the  observing  point 
0  by  the  Fresnel-Kirckoff  diffraction  formula  as  follows:  assume  the  occulter  S  lays 
on  a  plane  P  perpendicular  to  the  line  of  sight,  at  a  distance  D  from  the  observer 
(and  infinitely  far  from  the  source  of  light  in  the  plane  wave  approximation),  then 


05(0)  =  N 


g(  vW^+WTb^-D) 


(1  +  COS0)  dXdY, 


(1.1) 


where  P  —  S'  is  the  plane  of  the  occulting  screen,  6  the  angular  distance  from  the 
center  of  the  diffractor,  and  X  and  Y  the  Cartesian  distance  from  the  center  of 
the  screen  S  and  the  point  perpendicular  to  the  line  of  sight  on  the  plane  P  (see 
Figure  1.4).  This  intensity  is  normalized  to  05(0)  =  1  away  from  the  objects  by 
setting  =  l/2\i. 
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Figure  1.4:  Schematic  representation  of  an  occultation.  The  occulter  S  is  modeled  as  a 
circular  disk  laying  on  the  plane  P.  The  occultation  is  observed  by  an  observer  laying  at 
the  origin  0. 


When  the  occulter  approaches  the  line  of  sight,  D  3>  \/X‘^  +  then 


0^(0)  =  l-Xj  dX  dY. 


(1.2) 


If  the  occulter  is  a  circular  screen  with  radius  p  we  can  conveniently  move 
to  polar  coordinates  (/?,  0)  centered  on  the  center  of  the  occulter,  and  thus  the 
optical  path  difference  is  —  2  R  r  cos0  and  Equation  1.2 

can  be  expressed  in  terms  of  Bessel  functions  as  follows: 


ap(0) 

Jo{x) 


1 

TT 


2^  exp  ( 


i\D  J 
cosx  sint  dt 


exp 


^0 


R  dR 


(1.3) 


where  Jo  the  Bessel  function  of  order  0. 

In  units  of  Fresnel  scale,  F  =  a/AZ1/2,  the  integral  above  can  be  expressed 
as: 

ap(0)  =  1  +  f  Jo  {nr R)  R  dR  (1.4) 

Jo 
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and  by  ^  Jn+i{x))  =  x^^^Jn{x)  and  the  Lommel  functions: 


Unin,  =Y^  (-1)  (-)  ^n+2fc(7r/iZ/), 

k=0 

we  get  that  outside  of  the  geometrical  shadow,  or  for  r  > 

77r( -U  rp'^ 

ap(0)  =  1  +  fTrexp - ^ - -  {U2{p,r)  +  iUi{p,r)) ,  (1.5) 

while  inside  the  geometrical  shadow  {r  <  p),  using  we  have: 

j'jr( r‘^  /9^i 

ap(0)  =  exp - ^ - pUo{p,r) -iUi{p,r)) .  (1.6) 

Finally  it  follows  that,  the  measured  intensity  of  a  star  at  wavelength  A  is  described 
by: 

/ 

+  P<P 

Ip{p)  =  <  l  +  +  ^I(P>P)  P>P  _  (17) 

-2f/i(p,?7)  sin|(p2  +  r/2) 

^  +2f/2(p,p)cosf(p^ +  7]^) 

These  equations  describe  a  pattern  around  the  center  of  a  point  source 
star,  characterized  by  an  alternation  of  bright  and  dark  fringes  centered  on  the 
KBO.  During  the  transit  of  the  KBO  along  the  line  of  sight  this  translates  into  a 
modulated  lightcurve  (see  Figure  1.5).  This  basic  model  is  further  complicated  by 
the  hnite  size  of  the  star,  the  possibly  non-circular  shape  of  the  occulter  and  non- 
monocromatic  observations.  Some  of  these  points  will  be  addressed  in  the  following 
chapters.  Note  that  this  description  predicts  that  the  flux  in  the  center  of  the 
diffraction  pattern  will  be  exactly  equal  to  1  for  a  circular  occulter  of  size  smaller 
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than  the  Fresnel  scale  occnlting  at  6  =  0  impact  parameter.  This  point  in  the 
diffraction  pattern  is  called  Poisson  spot  in  honor  of  Poisson,  who  predicted  it 


Figure  1.5:  Diffraction  pattern  produced  by  a  D  =  3  km  KBO  at  42  AU  and  theoret¬ 
ical  diffraction  lightcurves  (in  magnitude  variation)  produced  by  the  observation  of  this 
occultation  at  different  impact  parameters  (right).  The  impact  parameters  are  marked  by 
horizontal  lines  on  the  left  panel. 

The  presence  of  diffraction  effects  in  the  event  of  an  occnltation  is  welcome 
to  the  observer  for  two  reasons.  The  transit  of  an  onter  solar  system  object  along  the 
line  of  sight  is  a  very  brief  event.  The  relative  velocity  of  the  objects  is  dominated 
by  the  velocity  at  which  the  earth  orbits,  ve  ~  30  km  s“^,  and  depending  npon 
the  observing  angle  (the  elongation)  TNOs  wonld  transit  across  the  line  of  sight 
at  a  speed  of  a  few  to  ~  25  km/s.  Thns  a  TNO  of  diameter  D  =  1  km  wonld 
transit  in  front  of  a  backgronnd  star  in  0.04-1  sec.  This  is  a  fast  rate  for  precision 
astrophysical  observations,  at  the  limits  of  feasibility  for  gronnd  based  snrveys. 
Diffraction  however  assnres  that  the  physical  size  of  the  event  is  no  smaller  than  the 
Fresnel  scale,  which  at  the  closer  end  of  the  Kniper  belt  is  abont  1  km  in  visible 

light.  In  a  diffraction  dominated  occnltation,  the  overall  flnx  rednction  is  dominated 

^Ironically  the  prediction  served  Poisson,  a  strong  believer  in  the  particle  nature  of  light,  to 
confute  Fresnel’s  theory  of  light  as  waves,  since  at  the  time  the  brightening  at  the  center  of  a 
shadow  had  not  been  observed. 


13 


by  the  size  of  the  KBO,  while  the  duration  of  the  event  depends  upon  the  relative 
velocity  Urei  and  the  size  of  the  diffraction  pattern  H.  The  relative  velocity  of  KBO 
can  be  approximated  to: 


^rel 


(1.8) 


where  Ae  the  distance  of  the  Earth  from  the  Sun,  and  e  the  angle  from  opposi¬ 
tion  (Liang  et  ah,  2004;  Nihei  et  ah,  2007).  We  dehne  the  cross  section  of  the  event 
H  as  the  diameter  of  the  first  Airy  ring:  the  first  (and  brightest)  bright  fringe  in 
the  occultation  pattern  (Born  &  Wolf,  1980).  H  it  is  limited  by  the  Fresnel  scale  for 
sub-kilometer  KBOs  and  by  the  size  of  the  object  for  large  KBOs  as  follows  (Nihei 
et  ah,  2007): 


H 


(2 


+  Ae, 


(1.9) 


where  the  additive  term  A6  accounts  for  the  hnite  angular  size  of  the  star.  When 
observing  at  opposition,  the  relative  velocity  Urei  of  an  object  orbiting  the  Sun  at 
40  AU  is  about  25  km  and  the  typical  duration  of  an  occultation  by  sub-kilometer 
KBOs  is  ~  0.2  s. 

Furthermore  the  occultation  features  contain  information  about  the  system 
that  generated  the  event,  with  potential  for  disentangling  the  size,  the  distance  and 
shape  of  the  occulter  (although  much  of  this  information  is  not  yet  accessible  in 
present  occultation  surveys).  Figure  1.6  shows  the  region  of  diameter  and  semi¬ 
major  axis  space  where  most  TNOs  reside.  The  current  limit  of  direct  observations 
is  shown  at  i?  =  30.  The  Fresnel  scale  is  also  shown:  all  objects  that  lay  below 
the  dashed  line  would  generate  diffraction  phenomena  during  an  occultation.  The 
angular  size  for  a  magnitude  K  =  13  GOV  star  is  shown:  as  the  angular  size  of  the 
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Figure  1.6:  Diameter  versus  semi-major  axis.  The  black  diamonds  show  the  known  KBOs. 
The  triangles  indicate  several  of  the  larger,  well-known  outer  Solar  System  objects  at 
their  semi-major  axes.  The  solid  lines  indicate  contours  of  constant  brightness  in  reflected 
sunlight,  assuming  an  albedo  value  of  0.04.  The  long  dashed  line  shows  the  Fresnel  scale 
as  a  function  of  distance  assuming  A  =  650  nm.  Occultations  by  objects  below  this  line 
are  diffraction  dominated.  The  dotted  line  is  the  angular  size  of  a  V=13  GOV  star  as  a 
function  of  distance.  The  limit  of  direct  observations  is  shown  at  i?  =  30. 

star  becomes  comparable  to  the  Fresnel  scale  the  diffraction  features  are  smoothed 
out  (see  Chapter  3  and  Nihei  et  ah  2007). 


1.4  The  story  told  by  the  small  KBOs 

Probing  the  very  small  {D  <  10  km)  region  of  the  KBO  size  spectrum 
and  the  regions  of  the  solar  system  outside  of  the  Kuiper  belt  could  have  profound 
consequences  on  our  understanding  of  the  formation  and  on  the  evolution  of  the 
solar  system. 

The  Kuiper  belt  has  been  shaped  by  accretion  and  disruption  processes 
throughout  the  history  or  the  solar  system.  With  small  orbital  eccentricities  the 


15 


relative  velocities  of  the  objects  in  the  early  Kuiper  belt  were  sufficiently  low  to 
allow  accretion  processes  to  form  kilometer  and  much  larger  objects.  Later,  when 
the  velocity  dispersion  increased  as  the  KBO  population  was  stirred  up  by  the  grav¬ 
itational  effects  of  the  larger  planets  and  planetoids,  only  large  objects  were  able  to 
continue  growing  through  impacts,  whereas  collisions  among  smaller  bodies  resulted 
in  disruption.  The  details  of  these  processes  depend  on  the  internal  strength  of  the 
KBOs  and  on  the  orbital  and  dynamical  evolution  of  the  gas  giant  planets.  The 
size  distribution  of  KBOs,  therefore,  contains  information  on  the  internal  structure 
and  composition  of  the  KBOs  -  and  hence  information  on  the  location  and  epoch 
in  which  they  formed  -  and  about  planetary  migration  (Pan  &  Sari,  2005;  Kenyon 
&  Bromley,  2004;  Kenyon  et  ah,  2008,  and  references  therein). 

Direct  observations  have  detected  KBOs  as  faint  as  magnitude  R  ~  28.2 
(Bernstein  et  ah,  2004),  which  corresponds  to  about  24  km  in  diameter  assuming 
a  4%  albedo.  The  KBO  size  distribution  can  be  characterized  using  its  brightness 
distribution.  The  latter  is  well  described  by  a  power  law  S(<  R)  =  deg“^, 

with  an  index  a  =  0.6  and  i?o  =  23  (Fraser  &  Kavelaars  2008,  Fuentes  &  Holman 
2008)  for  objects  brighter  than  about  R  =  25,  or  D  ~  100  km.  This  is  the  re¬ 
gion  of  the  size  spectrum  which  reflects  the  early  history  of  agglomeration.  Kenyon 
&  Windhorst  (2001)  pointed  out  that  the  intensity  of  the  infrared  Zodiacal  back¬ 
ground  sets  limits  on  the  extrapolation  of  a  straight  power  law  to  smaller  sizes. 
The  relatively  shallow  size  distribution  of  Jupiter  Family  Comets  (JFCs,  Tancredi 
et  ah  2006),  which  are  believed  to  originate  in  the  Kuiper  belt,  and  the  cratering 
of  Triton  observed  by  Voyager  2  (Stern,  1996),  all  point  to  a  flatter  distribution 
for  small  KBOs®.  In  2004  evidence  surfaced  that  a  break  in  the  power  law  occurs 

®The  relationship  between  the  cratering  of  Triton  and  the  Kuiper  belt  size  distribution  is 
questioned  in  Schenk  &  Zahnle  (2007). 
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at  a  diameter  larger  than  10  km:  Bernstein  et  al.  (2004)  condncted  deep  Hnbble 
Space  Telescope  observations  with  the  Advanced  Camera  for  Surveys  which  led  to 
the  discovery  of  only  3  new  objects  fainter  than  R  =  26,  about  4%  of  the  number 
expected  from  a  single  power  law  distribution  extrapolated  to  i?  =  29.  While 
this  work  remains  the  state  of  the  art  for  deep  direct  surveys  of  the  outer  solar  sys¬ 
tem,  recent  campaigns  have  observed  many  more  faint  objects  down  to  magnitude 
R  =  27,  which  with  the  assumption  of  a  4%  albedo  corresponds  to  about  40  km  in 
diameter^  (Fraser  &  Kavelaars  2008,  Fuentes  &  Holman  2008,  Fuentes  et  al.  2009b, 
and  Fraser  &  Kavelaars  2009).  These  recent  data  allowed  them  to  locate  a  break  in 
the  power  law  size  distribution  in  the  diameter  range  H  =  30  —  120  km. 

The  range  of  the  size  spectrum  of  Kuiper  belt  objects  (KBOs)  between  tens 
of  kilometers  and  meters  in  diameter  is  particularly  interesting  as  models  predict 
here  the  occurrence  of  transitions  between  different  strength  and  gravitation  regimes 
that  would  leave  a  signature  in  the  size  distribution  (Pan  &  Sari  2005,  Kenyon  & 
Bromley  2004,  Benavidez  &  Campo  Bagatin  2009,  and  references  therein).  Occulta- 
tion  surveys  allow  us  to  reach  farther  then  the  current  limits  of  direct  observations, 
and  into  this  very  region  of  interest,  and  they  are  the  only  observational  method 
presently  expected  to  be  able  to  detect  such  small  objects  in  the  outer  solar  system. 

While  occultation  surveys  were  first  proposed  in  1976,  only  recently  have 
results  been  reported.  This  observational  technique  requires  sub-second  photomet¬ 
ric  measurements  which  have  only  recently  become  possible.  Chang  et  al.  (2006) 
conducted  a  search  for  KBO  occultations  in  the  archival  Rossi  X-ray  Timing  Ex¬ 
plorer  (RXTE)  observations  of  Scorpius-Xl,  the  brightest  X-ray  source  in  the  sky. 
RXTE  is  a  satellite  dedicated  to  the  observation  of  X-ray  astronomical  sources,  able 

^The  magnitude  of  KBOs  is  converted  into  diameter  by  assuming  a  nominal  4%  albedo  through¬ 
out  this  work,  except  for  the  Fraser  &  Kavelaars  2008  result  where  the  authors  assumed  an  albedo 
of  6%. 
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to  provides  high  cadence  (>  fisec)  time  series  of  X-ray  sources.  Chang  et  ah  (2006) 
explored  nearly  90  hours  of  Sco-Xl  data  collected  between  1996  to  2002  by  RXTE, 
and  reported  a  surprisingly  high  rate  of  occultation-like  phenomena:  dips  in  the 
lightcurves  compatible  with  occultations  by  objects  between  10  and  200  m  in  di¬ 
ameter.  Jones  et  al.  (2008)  showed  that  most  of  the  dips  in  the  Sco-Xl  lightcurves 
may  be  attributed  to  artihcial  effects  of  the  response  of  the  RXTE  photo-multiplier 
after  high  energy  events,  such  as  strong  cosmic  ray  showers.  Only  12  of  the  original 
58  candidates  cannot  be  ruled  out  as  artifacts,  but  are  hard  to  conhrm  as  events 
(Jones  et  ah,  2008;  Chang  et  ah,  2007;  Liu  et  ah,  2008).  New  RXTE/PCA  data  of 
Sco  X-1  provided  a  less  constraining  upper  limit  to  the  size  distribution  of  KBOs 
(Liu  et  ah,  2008). 

Several  groups  have  conducted  occultation  surveys  in  the  optical  regime. 
Roques  et  al.  (2006)  and  Bickerton  et  al.  (2008)  independently  observed  narrow 
helds  at  45  Hz  and  40  Hz,  respectively,  with  frame  transfer  cameras.  Such  cameras 
allowed  them  to  obtain  high  signal-to-noise  ratio  (SNR)  fast  photometry  on  two  stars 
simultaneously.  Both  surveys  expect  a  very  low  event  rate  due  to  the  limited  number 
of  stars  and  the  limited  exposure,  and  neither  survey  has  claimed  any  detection  of 
objects  in  the  Kuiper  belt  at  this  time®.  An  upper  limit  for  KBOs  with  D  >1  km 
was  derived  by  Bickerton  et  al.  (2008)  by  combining  the  non-detection  result  of  the 
surveys  of  Chang  et  al.  (2007),  Roques  et  al.  (2006),  and  Bickerton  et  al.  (2008). 

In  my  graduate  studies  I  participated  in  two  campaigns  to  detect  occul¬ 
tation  events  in  star  time-series.  This  effort  is  described  in  the  following  chapters. 
Chapters  2-4  describe  the  effort  on  the  determination  of  the  size  distribution  of 
KBOs.  The  survey  I  conducted  at  the  MMT  with  the  Megacam  imager  is  described 
in  Chapter  2.  TAOS  (Taiwanese  American  Occultation  Survey)  is  a  dedicated  auto- 
®Roques  et  al.  report  3  possible  occultations  from  objects  outside  of  the  Kuiper  belt. 
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mated  multi-telescope  survey  (Leliner  et  al.,  2009b).  TAOS  reported  no  detections 
but  placed  the  strongest  upper  limit  to  date  to  the  surface  density  of  small  KBOs, 
which  is  reported  in  Zhang  et  ah  (2008).  My  work  on  the  hrst  3.75  years  of  TAOS 
data,  a  substantially  larger  dataset  than  the  one  used  in  Zhang  et  al.  (2008),  is  de¬ 
scribed  in  Chapters  3  and  4.  Chapter  5  describes  a  search  in  progress  for  Sedna-like 
and  scattered  disk  objects  in  the  TAOS  data.  Finally  I  summarize  the  conclusions 
and  discuss  perspectives  for  the  future  of  this  work. 
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Chapter  2 


The  sub-km  end  of  the  Kuiper 
Belt  size  distribution 


We  conducted  a  search  for  occultations  of  bright  stars  by  Kuiper  belt  Ob¬ 
jects  (KBOs)  to  estimate  the  density  of  sub-km  KBOs  in  the  sky.  We  report  here  the 
hrst  results  of  this  occultation  survey  of  the  outer  solar  system  conducted  in  June 
2007  and  June/ July  2008  at  the  MMT  Observatory  using  Megacam,  the  large  MMT 
optical  imager.  We  used  Megacam  in  a  novel  shutterless  continuous-readout  mode 
to  achieve  high  precision  photometry  at  200  Hz,  which  with  point-spread  function 
convolution  results  in  an  effective  sampling  of  ~30  Hz.  We  present  an  analysis  of  220 
star  hours  at  signal-to-noise  ratio  of  25  or  greater.  The  survey  efficiency  is  greater 
than  10%  for  occultations  by  KBOs  of  diameter  D  >  0.7  km,  and  we  report  no  de¬ 
tections  in  our  dataset.  We  set  a  new  95%  conhdence  level  upper  limit  for  the  surface 
density  Ti]\f{D)  of  KBOs  larger  than  1  km:  Ti]\f{D  >  1  km)  <  2.0  x  10®  deg“^,  and 
for  KBOs  larger  than  0.7  km  >  0.7  km)  <  4.8  x  10®  deg“^.^ 

^This  work  appeared  in  The  Astronomical  Journal,  Volume  138,  Issue  2,  pp.  568-578  (2009), 
with  title  A  Search  for  Occultations  of  Bright  Stars  by  Small  Kuiper  Belt  Objects  Using  Megacam 
on  the  MMT  and  it  was  co-authored  by  P.  Protopapas  (Harvard-Smithsonian  Center  for  Astro¬ 
physics,  Initiative  in  Innovative  Computing  at  Harvard),  B.  A.  McLeod  (Harvard-Smithsonian 
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2.1  Introduction 


The  survey  I  report  here  was  conducted  using  Megacam  (McLeod  et  ah, 


2006,  Figure  2.2)  at  the  6.5  m  MMT  Observatory  at  Mount  Hopkins,  Arizona.  The 


use  of  Megacam  in  continuous-readout  mode  (see  Section  2.2)  on  a  held  of  view 


of  24'  X  24'  allowed  us  to  monitor  over  ~  100  stars  at  200  Hz  over  the  course  of 


two  observational  campaigns  conducted  in  June  2007  and  June-July  2008.  This 


peculiar  use  of  a  conventional  CCD  camera  allowed  us  to  reach  the  high  speed 


photometric  sampling  necessary  to  detect  occupations  by  small  outer  solar  system 


objects  without  compromising  the  number  of  star  targets  monitored.  Our  survey  is 
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Figure  2.1:  A  simulated  diffraction  pattern  (left  panel)  generated  by  a  spherical  D  =  1km 
KBO  occulting  a  magnitude  12  FOV  star.  The  MMT/Megacam  system  bandpass  (Sloan 
r’  filter  and  camera  quantum  efficiency)  is  assumed.  The  size  of  the  KBO  and  the  size  of 
the  Airy  ring  -  a  measure  of  the  cross  section  of  the  event  -  are  shown  for  comparison. 
The  right  panel  shows  the  diffraction  signature  of  the  event  (assuming  central  crossing: 
impact  parameter  b  =  0)  as  a  function  of  the  distance  to  the  point  of  closest  approach 
(bottom  scale).  The  top  scale  shows  the  time-line  of  the  event  assuming  an  observation 
conducted  at  opposition  (relative  velocity  Vrei  =  25  km  s“^).  The  occupation  is  sampled 
at  200  Hz  (dashed  line),  and  at  30  Hz,  the  effective  sampling  rate  after  taking  PSF  effects 
into  account  (solid  line,  see  Section  2.3). 

Center  for  Astrophysics),  C.  R.  Alcock  (Harvard-Smithsonian  Center  for  Astrophysics),  M.  J.  Hol¬ 
man  (Harvard-Smithsonian  Center  for  Astrophysics),  M.  J.  Lehner  (Institute  of  Astronomy  and 
Astrophysics  -  Academia  Sinica,  University  of  Pennsylvania,  Harvard-Smithsonian  Center  for  As¬ 
trophysics) 
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sensitive  to  occultations  by  outer  solar  system  objects  of  diameter  D  ~  700  m  or 
larger. 

In  our  survey,  the  bandpass  of  the  observation  is  centered  near  A  =  500  nm 
and,  at  distance  A  ~  40  AU,  the  Fresnel  scale  F  is  F  ~  1.2  km.  Any  occultation 
caused  by  objects  in  the  Kuiper  belt  of  a  few  kilometers  in  diameter  or  smaller 
will  exhibit  prominent  diffraction  effects.  Figure  2.1  shows  a  diffraction  pattern 
generated  by  a  H  =  1  km  KBO  (left)  and  the  diffraction  feature  that  would  be 
imprinted  in  a  star  lightcurve  observed  by  our  system  (right). 

We  report  no  detections  in  220  star  hours.  Our  MMT  survey  is  designed 


Figure  2.2:  Megacam  focal  plane  (McLeod  et  ah,  2006).  A  thick  rectangle  outlines  a  single 
CCD  in  the  9x4  CCDs  mosaic.  Two  halves  of  each  CCD  (thin  rectangles)  are  read  into  two 
separate  amplifiers;  each  amplifier  generates  a  separate  output  image  in  our  observational 
mode.  The  x  and  y  axis,  as  they  would  appear  in  a  resulting  image,  are  also  shown. 
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to  be  complementary  to  TAOS  and  to  reach  smaller  size  limits,  and  unlike  TAOS 
it  would  allow  us  to  estimate  the  size  of  a  detected  occulting  KBO.  We  expect 
further  work  on  adaptive  photometry  and  de-trending  to  signihcantly  improve  our 
sensitivity,  perhaps  allowing  us  to  detect  KBOs  as  small  a.s  D  >  300  m.  I  discuss 
the  improvements  we  are  developing  on  this  analysis  in  Section  6.1.  The  preliminary 
analysis  we  present  here  allows  us  to  derive  upper  limits  for  objects  D  >  700  m. 

In  the  next  section  I  describe  the  novel  observational  mode  adopted  for 
this  survey.  In  Section  2.3  I  describe  the  data  acquired  and  analyzed  for  this  paper. 
Details  of  the  data  extraction  and  reduction,  which  required  custom  packages,  are 
addressed  in  the  same  section.  Section  2.4  describes  the  characteristics  of  the  noise 
of  our  current  datasets,  and  our  noise  mitigation  approach.  Section  2.5  describes 
the  detection  algorithm.  In  Section  2.6  I  derive  our  upper  limit  to  the  density  of 
KBOs.  I  also  compare  in  detail  the  achievements  of  our  survey  to  those  of  previous 
surveys. 


2.2  Fast  Photometry  with  a  Large  telescope:  The 
Continuous  —  Readout  Mode 

Achieving  sub-second  photometric  sampling  is  a  challenge  in  optical  as¬ 
tronomy.  CCD  cameras  can  perform  fast  photometric  observations  by  reading  out 
small  sub-images,  limiting  the  observations  to  very  small  portions  of  the  sky  (e.g.. 
Marsh  &  Dhillon  2006,  Bickerton  et  ah  2008).  This  is  the  approach  adopted  by 
Roques  et  ah  (2006),  and  Bickerton  et  al.  (2008),  who  observed  two  stars  at  one 
time.  Due  to  the  rarity  of  occultation  events,  however,  one  would  want  to  maximize 
the  number  of  targets  and  the  total  exposure  to  increase  the  number  of  detections. 
TAOS  achieves  sub-second  photometric  observation  on  up  to  500  targets  with  the 
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Figure  2.3:  Conventional  stare  mode  image  (one  half  of  a  CCD)  of  one  of  our  fields  (bottom 
panel).  A  series  of  rows  from  continuous-readout  mode  (center  panel)  from  the  same  CCD 
and  field,  where  the  rows  are  stacked  together  in  a  single  image.  The  flux  profile  of  the 
central  row  of  this  segment  of  continuous-readout  data  is  plotted  in  the  top  panel. 

zipper  mode  readout  technique  (Lehner  et  ah,  2009b),  but  they  sample  at  <  5  Hz 
rate.  Our  continuous-readout  technique  allows  us  to  observe  the  entire  held  of  view 
of  the  camera  at  200  Hz. 

Megacam,  the  MMT  optical  imager,  is  a  mosaic  camera  comprising  thirty- 
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six  CCDs  -  each  with  an  array  of  2048  x  4608  pixels  -  with  a  24'  x  24'  field  of 
view  (Figure  2.2).  The  standard  readout  speed  of  each  CCD  is  0.005  sec/row  with 
2x2  binning.  For  this  survey,  we  operated  the  camera  in  shutterless  continuous- 
readout  mode;  that  is,  we  kept  the  shutter  open  while  scrolling  and  reading  the 
charges  at  the  standard  readout  speed,  tracking  the  sky  at  the  sidereal  rate.  Each 
star  is  represented  in  each  row  that  is  read  out  of  the  camera,  and  the  flux  from  a 
star  in  a  row  represents  a  photometric  measurement  of  that  star  sampled  at  200  Hz. 
Stacking  each  read  row  into  a  single  image  each  star  time-series  forms  a  streak  along 
the  readout  axis  (y— axis).  A  small  portion  of  our  data  is  shown  in  Figure  2.3. 

In  this  observational  mode  the  flux  from  the  sky  background  is  added 
continuously  as  the  charge  is  transferred  from  one  end  of  the  CCD  to  the  other,  so 
the  sky  is  exposed  2304  x  0.005  =  11.52  sec  for  every  0.005  sec  integration  on  each 
star  image  (where  2304  is  the  effective  number  of  rows  in  each  2x2  binned  CCD). 
In  this  mode  the  photon  limited  SNR  is  typically  ~  180  for  an  r'  magnitude  10  star. 

When  observing  multiple  targets  simultaneously  one  can  notice  that  the 
lightcurves  are  affected  by  common  fluctuations,  or  trends,  due  for  example  to 
weather  patterns  (Kim  et  ah,  2008,  and  references  therein).  In  our  observational 
mode,  however,  additional  flux  variations  are  caused  by  wind-induced  resonant  oscil¬ 
lations  of  the  telescope.  While  the  image  motion  along  the  x  axis  of  the  focal  plane 
(transverse  to  the  readout  direction)  can  be  resolved  (see  Section  2.3.1),  the  image 
motion  parallel  to  the  direction  of  the  CCD  readout  induces  an  effective  variation 
in  the  exposure  time  of  a  star  for  a  given  row.  These  fluctuations  are  common  to 
all  stars  in  the  field  (with  possible  position  dependencies)  and  therefore,  in  prin¬ 
ciple,  they  are  completely  removable.  We  discuss  the  de-trending  of  our  data  in 
Section  2.3.1.  Other  sources  of  noise  that  affect  continuous-readout  mode  data  are 
discussed  in  Section  2.4.  The  typical  duration  of  a  set  of  contiguous  data  was  10-15 
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Figure  2.4:  Top:  power  spectrum  of  one  of  our  lightcurves  before  and  after  de-trending  the 
lightcurve  to  remove  noise  (see  Section  2.3.1).  Bottom:  power  spectrum  of  the  occultation 
time-series  for  a  1  km  KBO  at  40  AU  occulting  a  FOV  V  =  12  star. 

minutes  (after  which  the  data  load  on  the  buffer  would  become  prohibitive).  For 
each  ampliher,  a  single  FITS^  hie  is  created  wherein  all  of  the  rows  read  out  during 
a  data  run  are  stored  as  a  single  image.  For  a  typical  run  each  FITS  output  image 
contains  lOOK  to  130K  rows,  corresponding  to  about  150-200  Mb  of  data. 

2.3  Data 

We  selected  observing  helds  within  2.8°  of  the  ecliptic  plane,  where  the 
concentration  of  KBOs  is  highest  (Brown,  2001).  In  order  to  maximize  the  number 
of  targets  we  selected  our  helds  at  the  intersection  of  the  ecliptic  and  galactic  planes 
(RA  ~  lO'^OOOO®,  Dec  ~  — 21°00'00").  We  conducted  our  observations  in  June-July, 
when  our  helds  were  near  opposition  (elongation  angle  e  =  180°)  and  the  relative 
velocity  of  the  KBOs  is  highest  (Roques  et  al.,  1987;  Nihei  et  al.,  2007;  Bickerton 
^Flexible  Image  Transport  System,  http://fits.gsfc.nasa.gov/. 
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Table  2.1.  MMT/Megacam  survey  observed  fields 


RA 

Dec 

A'’ 

(deg) 

e  range’’ 
(deg) 

17’’0000" 

-21°15'00" 

1.5 

174-160 

17’’1500" 

-20°15'00" 

2.8 

176-163 

18'’0000" 

-21°15'00" 

2.2 

171-173 

18'’0000" 

-21°30'00" 

1.9 

171-173 

18'’0000" 

-21°45'00" 

1.7 

172-173 

19'’0000" 

-22°00'00" 

0.7 

158-172 

’^ecliptic  longitude 
’’range  of  elongation  angles 

et  al.,  2009),  thus  maximizing  the  event  rate  per  target  star.  Pointing  information 
for  our  fields  is  summarized  in  Table  2.1.  The  RA  and  Dec  of  each  observed  field 
are  listed  together  with  the  ecliptic  latitude  (A)  and  a  maximum  range  of  elongation 
angles  at  which  the  filed  might  have  been  observed. 

We  also  observed  control  fields.  These  were  chosen  on  the  galactic  plane 
at  a  high  ecliptic  latitude;  we  expect  a  negligible  rate  of  occultations  by  KBOs  in 
these  fields.  These  data  allow  us  to  assess  our  false  positive  rate.  Since  we  report  no 
detections  the  analysis  of  these  fields  is  not  discussed  further  in  this  paper.  All  of 
our  observations  were  conducted  in  Sloan  r'  filter  (Smith  et  ah,  2002).  A  set  of  about 
7  hours  on  target  fields  was  collected  in  5  half  nights  in  June  2007  and  a  similar 
number  of  hours  was  collected  on  control  fields.  A  set  of  about  7  hours  on  target 
fields  and  about  6  hours  on  control  fields  was  collected  in  7  half  nights  in  June-July 
2008.  Out  of  the  2007  dataset  100.61  star  hours  at  SNR  >  25  are  considered  in 
this  paper.  From  the  2008  dataset  we  use  here  118.93  star  hours.  Information  on 
our  dataset  is  summarized  in  Table  2.2.  The  minimum  signal-to-noise  ratio  of  25,  is 
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Table  2.2.  MMT/Megacam  survey  data  set  parameters. 


Start  Date 

End  Date 

Exposure  at  SN>  25 

Number  of  lightcurves  with  SN>  25 
Number  of  Photometric  Measurements 

2007  June  6 

2007  June  10 
100.61  star-hours 
990 

7.2  X  10^ 

Start  Date 

End  Date 

Exposure  at  SNR>  25 

Number  of  lightcurves  with  SNR>  25 
Number  of  Photometric  Measurements 

2008  June  27 

2008  July  1 

118.93  star-hours 
527 

8.5  X  10^ 

chosen  arbitrarily:  25  is  the  minimum  SNR  of  the  surveys  of  Roques  et  al.  (2006) 
and  Bickerton  et  al.  (2008).^  A  SNR  25  limits  our  sensitivity  to  fluctuations  greater 
than  4%.  An  occultation  of  a  magnitude  12  FOV  star  by  a  KBO  oi  D  =  400  m 
diameter  would  produce  a  4%  effect.  Our  efficiency  tests,  however,  revealed  our 
sensitivity  rapidly  drops  below  10%  for  objects  smaller  than  D  =  700  m,  due  to 
residual  non-Gaussianity  in  our  time-series  photometric  data.  We  discuss  this  in 
Section  2.4. 


2.3.1  Data  extraction  and  reduction 

Extraction 

Custom  algorithms  have  been  developed  for  the  data  extraction  and  re¬ 
duction.  For  each  field  a  preliminary  stare  mode  (conventional)  image  is  collected 
before  each  series  of  high-speed  runs.  At  the  beginning  of  our  analysis  the  stare 

^Note  however  that  this  SNR  level  is  obtained  here  for  200  Hz,  whereas  Roques  et  al.  (2006) 
and  Bickerton  et  al.  (2008)  observed  at  45  Hz  and  40  Hz. 


mode  image  is  analyzed  using  SExtractor  (Bertin  &  Arnouts,  1996)  to  generate  a 
catalog  of  bright  sources.  This  catalog  is  used  to  identify  the  initial  position  and 
brightness  of  each  star  in  the  focal  pane.  In  order  to  analyze  the  continuous  readout 
data,  we  hrst  determine  the  sky  background  for  each  CCD  and  each  row.  To  do 
so  we  calculate  the  mean  of  the  flux  counts  in  each  row  after  removing  the  mea¬ 
surements  that  are  three  a’s  or  more  above  the  mean  (Sa-clipping)  iteratively  until 
the  mean  converges.  This  removes  most  of  the  pixels  in  the  row  containing  flux 
from  resolved  stars.  Next,  a  subset  of  stars  that  are  bright  and  isolated  is  selected 
from  the  stare-mode  catalog  and  used  to  determine  the  ^-displacement  of  the  focal 
plane.  The  focal  plane  is  split  into  two  halves,  9x2  chips  each,  that  are  analyzed 
separately.  We  select  eight  stars,  two  near  each  of  the  four  corners  of  each  half-focal 
plane.  This  allows  us  to  characterize  the  global  motion  of  the  targets  even  in  the 
presence  of  small  rotational  modes  or  spatial  dependency  (see  Section  2.4).  For 
each  star  (x),  and  at  each  time-stamp  (t),  we  calculate  /i*(t)  and  a*(t),  respectively 
the  centroid  offset  from  the  original  position  and  the  standard  deviation  of  the  star 
image,  assuming  a  Gaussian  prohle.  Note  that,  for  a  given  time-stamp,  flux  from 
different  stars  will  appear  on  different  rows  due  to  the  y-positions  of  the  stars  on 
the  focal  plane.  A  1-D  Gaussian 


F, 


h  exp 


V  2a2(t)  ; 


+  hg 


(2.1) 


(where  F*  is  the  total  star  flux,  /*  the  flux  at  the  peak  and  Ag  the  sky)  is  £t  for 
each  of  the  eight  stars  to  each  row  of  the  star-streak.  Thus  the  ^-displacement  Ji{t) 
for  all  the  stars  in  the  held  at  time-stamp  t  is  estimated  to  be  the  weighted  average 
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of  the  star  displacements: 


m  =  ^ - ,  (2.2) 

E 
★  =  1 


where  p*(to)  is  the  star  initial  x-position  and  cn*  is  the  weight  used  for  that  star. 

In  order  to  weight  our  average  we  use  the  correlation  of  the  entire  x- 
displacement  time-series  p*  with  respect  to  the  rest  of  the  star  set: 


(2.3) 

(2.4) 


where  is  the  variance  of  the  displacement  throughout  the  duration  T  of  the  time- 
series.  The  weight  cu*  is  the  square  of  the  Pearson’s  correlation  coefficient  (Rice, 
2006,  pag.  406),  a  measure  of  the  correlation  of  the  displacement  time-series  for 
one  star  with  the  other  seven.  All  star  lightcurves  in  the  field  are  then  extracted 
by  aperture  photometry  adjusting  time-stamp  by  time-stamp  the  center  of  the 
aperture  according  to  the  x-motion  derived  in  this  stage,  and  with  a  hxed  aperture 
size  which  is  proportional  to  the  average  FWHM  in  the  run.^ 


De-trending 

The  lightcurves  thus  extracted  show  evident  semi-periodic  and  quasi-  sinu¬ 
soidal  flux  variations  that  can  be  associated  with  oscillatory  modes  of  the  telescope 

^We  attempted  to  extract  the  lightcurves  with  both  fixed  aperture  size  and  variable  aperture 
size,  using  the  FWHM  calculated  by  Gaussian  fitting  as  a  point  by  point  estimator  of  the  aperture 
size.  The  fixed  aperture  extraction  proved  to  be  more  reliable  than  the  variable  aperture  extraction, 
which  induced  further  noise  in  our  lightcurves. 
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Figure  2.5:  Image  motion  and  PSF  over  time:  mean  of  the  x  displacements  for  eight 
bright  isolated  stars,  at  the  four  corners  of  the  half-focal  plane  for  two  data  runs  (top  left 
and  right  panels).  PSF  width  from  the  Gaussian  fit  averaged  over  the  same  set  of  stars 
(bottom  left  and  right).  On  the  left  we  used  the  same  run  used  to  generate  Figure  2.8. 
The  arrow  points  to  the  displacement  feature  marked  in  Figure  2.8.  On  the  right  the 
x-displacement  and  the  PSF  width  for  another  run,  with  the  first  0.5  seconds  shown  on 
the  left  at  higher  time  resolution.  Note  how  in  the  second  run  the  x-displacements  are 
less  prominent  (note  the  different  y  scale)  but  the  amplitude  of  the  variability  of  the  PSF 
is  larger. 

in  the  y  direction.  In  particular,  a  Fourier  analysis  generally  reveals  two  strong 
modes,  roughly  consistent  among  runs,  one  with  period  near  0.04  seconds  and  the 
other  near  0.5  seconds.  Fourier  spectra  for  one  of  our  lightcurves,  before  and  after 
processing  it,  are  shown  in  Figure  2.4  (top).  Because  these  fluctuations  affect  the 
whole  CCD  plane,  they  are  common  to  all  stars  and  can  be  removed  to  achieve 
greater  photometric  precision.  We  now  want  to  identify  and  remove  these  trends 
from  our  lightcurves,  a  process  that  we  call  de-trending. 

The  general  algorithm  we  used  for  de-trending  is  described  in  Kim  et  ah 
(2008).  The  method  takes  advantage  of  the  correlation  among  lightcurves  to  extract 
and  remove  common  features.  Since  we  can  identify  distinct  semi-periodic  modes 
we  de-trend  high  and  low  frequencies  separately  (typically  u  >  lOHz  and  u  <  lOHz). 

We  hrst  smooth  the  lightcurves,  to  remove  all  but  the  frequencies  that  we 
want  to  de-trend,  by  applying  a  low-pass  or  high-pass  hlter.  We  then  select  a  subset 
of  Nr  template  lightcurves  (/,-)  that  show  the  highest  correlation  in  the  lightcurve 
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features.  is  typically  about  15.  A  master  trend  lightcurve  r  is  generated  as  the 
weighted  average  of  the  normalized  template  lightcurves: 

E  /rj(t)/(/r,i) 

"W  =  F 

E 

i=i 

where  the  notation  (/tj)  denotes  the  mean  flux  of  frjit)  over  the  duration  T  of  the 
lightcurve,  and  the  weight  is  the  variance  of  the  lightcurve  in  time;  r(t)  has 

mean  value  of  unity  and  it  represents  the  correlated  fluctuations  in  all  lightcurves. 

The  main  trend  is  physically  associated  with  an  over-under  exposure  phe¬ 
nomenon  due  to  global  image  motion  along  the  y  axis,  which  causes  the  effective 
exposure  time  to  vary  (see  Section  2.2),  therefore  scaling  the  flux.  In  order  to  remove 
these  common  trends  we  divide  point  by  point  the  flux  of  each  original  lightcurve  / 
by  the  trend  master  lightcurve.  To  improve  the  de-trending  effectiveness  we  allow 
a  free  multiplicative  factor  A/  (a  scaling  factor)  for  each  lightcurve  as  follows: 

U.A,W  =  f(t)  +  1  ;  (2,6) 

fd,Af  is  the  de-trended  lightcurve. 

We  optimize  our  de-trending  by  selecting  Af  to  minimize  the  variance  of 
the  de-trended  lightcurve  /d  with  respect  to  /c  =  /  —  /s  +  (/s) ,  which  is  the  original 
lightcurve  cleaned  of  the  frequency  to  be  de-trended.  We  apply  a  high-pass  (low- 
pass)  hlter  to  /  to  obtain  if  we  want  to  de-trend  the  low  (high)  frequencies.  Af 
is  then  optimized  by  setting: 

^E(/d,A,(«)  -  (/c))"  =  0.  (2.7) 

f  t=l 
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which  minimizes  the  second  moment  of  the  de-trended  lightcurve  with  respect  to  fc- 
The  optimal  value  of  A f  can  be  calculated  analytically. 

We  set  no  constraints  on  Af,  and  for  all  of  our  runs  the  optimal  values 
of  A f  proved  to  be  close  to  1  (which  is  what  we  expect  in  the  presence  of  global 
trends)  except  for  pathological  cases  where  the  flux  of  the  star  was  buried  in  noise 
and  the  raw  and  de-trended  SNR  were  extremely  low.  These  lightcurves  would  not 
pass  SNR  cuts  and  were  never  considered  in  any  of  our  analysis. 

Examples  of  the  results  obtained  by  our  de-trending  algorithm  are  dis¬ 
played  in  Figures  2.6  and  2.7.  In  Figure  2.6  the  top  two  panels  show  lightcurves 
for  two  independent  sources  in  our  field,  and  the  bottom  two  panels  show  the  same 
lightcurves  after  de-trending.  Note  that  the  top  star  is  ~  2.5  magnitudes  brighter 
than  the  other  and  this  is  reflected  in  the  lower  SNR  of  the  fainter  source  (bottom 
panel).  Figure  2.7  shows  one  of  our  lightcurve  before  (top)  and  after  de-trending 
(bottom).  The  raw  lightcurve  is  implanted  with  an  occultation  by  a  D  =  1  km 
KBO  occulting  a  V  =  9  FOV  star.  The  diffraction  feature  is  completely  lost  in  the 
trends  and  becomes  evident  only  after  de-trending.  In  the  bottom  panel  we  show 
the  lightcurve  de-trended  without  allowing  for  the  optimization  factor  Af  at  the  top 
(plotted  at  the  top  at  an  arbitrary  offset)  and  with  optimization  factor  Af  =  1.15 
for  the  low  frequencies  and  Af  =  1.05  for  the  high  frequencies,  shown  at  the 

bottom.  The  introduction  of  an  optimization  factor  improves  the  SNR  of  the  de¬ 
trended  lightcurve  from  SNR  =  30.0  to  SNR  =  30.7.  For  this  particular  run 
improvements  of  up  to  7%  in  SNR  were  achieved  by  optimizing  the  de-trending. 

Note  that,  while  we  used  smoothed  versions  of  our  lightcurves  to  identify 
the  trends  and  to  optimize  the  de-trending,  we  do  not  smooth  or  hlter  our  lightcurves 
to  improve  the  SNR,  thus  preserving  all  intrinsic  features  (including  potential  oc- 
cultations).  Figure  2.4  shows  the  power  spectrum  of  one  of  our  lightcurves  before 


33 


0  O.S  0.4  18  3  4 

t  (seconds) 


Figure  2.6:  Lightcurves  of  two  independent  stars  in  one  of  our  fields.  The  left-hand  plots 
show  a  0.5-second  chunk  of  the  time  series;  the  following  4-seconds  are  shown  on  the  right 
at  a  lower  time  resolution.  The  top  two  panels  show  the  lightcurves  before  de-trending. 
Common  modes  are  visible  at  multiple  time  scales.  The  bottom  two  panels  show  the 
lightcurves  after  de-trending.  The  top  lightcurve  is  the  same  used  in  Figure  2.4 

and  after  de-trending  it  (top).  The  power  spectrum  of  an  occultation  time-series 
generated  by  a  1  km  KBO  occulting  a  FOV  star  of  magnitude  V  =  12  is  shown 
in  the  bottom  panel.  Our  de-trending  greatly  reduced  the  power  at  all  frequencies: 
the  cumulative  power  for  this  particular  lightcurve  at  frequencies  z/  <  40  Hz  is  sup¬ 
pressed  by  a  factor  of  40.  Because  the  oscillations  are  not  perfectly  correlated  among 
our  stars  (see  Section  2.4)  some  residual  power  is  visible.  Smoothing  however  would 
would  signihcantly  reduce  the  strength  of  the  occultation  features,  that  show  power 
at  all  frequencies  z/  <  20  Hz. 
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Figure  2.7:  Raw  lightcurve  on  which  the  occultation  signature  of  a  1  km  KBO  occult¬ 
ing  a  magnitude  F  =  9  FOV  star  has  been  implanted  (top)  and  the  same  lightcurve 
after  de-trending  (bottom).  In  the  bottom  panel  the  top  lightcurve  is  de-trended  without 
optimization  (plotted  at  the  top  at  an  arbitrary  offset)  and  the  bottom  lightcurve  is  de¬ 
trended  with  optimization  factor  Af  =  1.15  for  the  low  frequencies  and  1.05  for  the  high 
frequencies. 

2.4  Residual  noise  in  the  time-series 

With  a  SNR  >  25  we  can  detect  fluctuations  of  a  few-percent.  In  an 
0.005  sec  exposure  the  flux  for  a  magnitude  r'  =  14  star  observed  by  Megacam  is 
about  10^  e~ ,  which  after  taking  into  account  the  contribution  to  the  noise  of  the 
background  should  lead  to  a  Poisson  limited  SNR  of  about  25.  While  we  were  able 
to  remove  a  large  portion  of  the  noise  that  originally  affected  our  data,  we  typically 
cannot  reach  the  Poisson-limit.  We  have  identihed  hve  possible  sources  of  noise  in 
our  data: 

•  Contamination  by  nearby  sources.  Overlap  of  stars  along  the  x  axis  (per¬ 
pendicular  to  the  read-out  direction)  within  a  chip,  causes  reciprocal  contam¬ 
ination  in  our  readout  mode  and  some  of  the  stars  in  our  helds  are  therefore 
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compromised  and  excluded  from  our  analysis.  Furthermore,  oscillations  of  the 
images  along  the  x  axis  causes  the  relative  distance  between  the  star-streaks 
to  change,  which  causes  occasional  merging.  Note  that  while  these  oscillations 
are  simultaneous  in  time  domain,  they  do  not  occur  in  the  same  row  in  the 
recorded  image.  In  each  row  the  star  images  of  two  objects  that  are  at  a  dif¬ 
ferent  y  position  on  the  CCD  plane  will  not  belong  to  the  same  time-stamp, 
therefore  the  oscillations  -  while  simultaneous  -  will  show  a  y  offset.  This  is 
shown  in  Figure  2.8.  The  merging  of  streaks  causes  artihcially  high  counts. 
Aperture  photometry  with  a  fixed  aperture  does  not  address  this  issue  properly 
and  fitting  photometry  on  individual  streaks  is  a  computationally  expensive, 
inefficient  method  which  is  also  unstable  in  the  presence  of  multiple  sources 
close  to  each  other. 

•  Unresolved  sources.  Sources  that  are  too  faint  to  be  visible  in  our  0.005 
sec  exposures  generate  a  diffuse  background.  For  the  data  in  Figure  2.3  the 
sky  level  calculated  as  the  Sa— clipped  mean  of  the  row  counts  is  140.5  ADUs. 
The  stare  mode  image  sky  level  was  48  ADUs  for  a  5  sec  exposure,  which 
would  lead  to  a  prediction  of  110  ADUs  for  our  0.005  x  2304  sec  effective 
continuous-readout  exposure.  The  discrepancy  is  due  to  the  presence  of  un¬ 
resolved  streaks  associated  with  faint  stars  across  the  held.  Summing  all  the 
counts  in  the  stare  mode  image  and  rescaling  by  the  exposure  time  of  each  row 
we  get  a  number  very  close  to  the  sum  of  all  counts  in  a  row  of  continuous- 
readout  data.  This  contamination  introduces  extra  Poisson  noise,  but  more 
importantly  it  introduces  non-Poissonian  noise  as  well,  since  the  unresolved 
sources  are  affected  by  the  same  trends  the  bright  stars  display.  Our  data 
shows  evidence  of  off-phase  correlation  that  might  be  induced  by  unresolved 
sources. 
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•  Positional  dependency  in  the  motion  and  trends.  While  we  treat  all  of 
the  stars  in  the  held  as  an  ensemble  that  moves  in  a  solid  fashion  along  the  x 
and  y  axes,  the  image  motion  might  also  have  a  rotational  component.  This 
would  lead  to  position  dependencies  in  the  motion  that  are  not  accounted  for 
by  our  aperture  centering  algorithm.  We  have  not  seen  evidence  of  dependency 
on  the  distance  to  the  center  of  the  focal  plane  in  either  motion  or  trends,  but 
we  cannot  exclude  that  occasional  rotational  modes  of  the  telescope  would 
occur.  Differential  image  motion  and  hux  huctuations  might  also  be  induced 
by  atmospheric  seeing.  Both  of  these  effects  might  cause  the  star-streak  to 
move  out  of  the  photometric  aperture  leading  to  artihcially  low  counts.  The 
aperture  size  must  be  chosen  to  be  such  that  errors  due  to  contamination  by 
nearby  sources  and  errors  due  to  streaks  exiting  the  aperture  are  simultane¬ 
ously  minimized.  Furthermore  in  the  presence  of  thin  clouds,  variations  in 
the  transparency  might  generate  trends  that  would  affect  different  sources  at 
different  times  as  the  clouds  move  across  the  image.  Positional  dependencies 
or  variations  in  transparency  might  contribute  to  the  off-phase  correlation  of 
our  data. 

•  Scintillation.  Young’s  scaling  law  (Young,  1967), 

0.09A^/^(sec  ZY-'^^exp{—^ 

^  ~  V2  At 

describes  the  error  due  to  the  low-frequency  component  of  scintillation,  with 
a  =  {AI/I)  and  where  A  is  the  telescope  aperture,  Z  the  angle  from  zenith,  h 
is  the  height  of  the  turbulence  layer,  and  AT  the  integration  time.  Competing 
effects  are  in  place  in  our  survey:  the  large  aperture  mitigating  the  noise, 
and  the  low  air  mass  contributing  to  noise  degradation.  Using  the  above 
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equation  and  representative  data  from  La  Palma  (Dravins  et  al.,  1998)  we 
estimate  that  the  noise  contribution  from  scintillation  is  aj  <  0.01,  i.e.,  not 
the  dominant  source  of  residual  noise.  Note  however  that  this  relation  holds 
for  integration  time  on  scales  of  seconds  or  longer.  When  including  the  effects 
of  high-frequency  scintillation  the  dependency  on  the  aperture  is  expected  to 
be  steeper: 

poo 

aj  (X  A~'^^^{sec  Z)^  /  C‘^{h)h‘^dh,  (2.9) 

Jo 

where  Cn  is  the  refraction  coefficient  for  the  turbulent  layer  (see  Dravins  et  al. 
1997  and  references  therein). 

As  compared  to  the  other  occnltation  surveys  the  term  associated  to  the  tele¬ 
scope  aperture  (A”’^/®)  in  the  SNR  variance  is  a  factor  20  lower  than  the  same 
factor  for  the  TAOS  survey,  4.5  times  lower  than  the  same  factor  for  Bickerton 
et  al.  (2008)  and  1.4  than  for  the  Roques  et  al.  (2006)  survey. 

•  Convolution  of  the  time  series  with  finite  PSF.  The  hnite  size  of  the  PSF 

(typically  two  to  three  pixels,  although  it  occasionally  was  as  large  as  seven) 
causes  consecutive  measurements  to  be  correlated.  This  effect  is  not  a  source 
of  noise  per  se,  but  it  changes  the  spectral  characteristics  of  the  noise.  The 
scale  of  this  phenomenon  shows  up  in  an  anto-correlation  analysis  with  high 
power  at  a  lag  of  abont  seven  pixels.  This  is  effectively  a  kernel  convolntion  of 
onr  time  series  that  smooths  the  signal,  including  possible  occnltation  signals, 
so  that  while  we  sample  the  images  at  200  Hz  we  would  expect  an  occnltation 
signatnre  to  be  effectively  sampled  at  ~  30  Hz  (see  Section  2.5).  Note  that  this 
is  close  to,  bnt  slightly  short  of,  the  critical  Nyqnist  sampling  for  occnltations 
dominated  by  diffraction  (Bickerton  et  ah,  2009). 
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Figure  2.8:  A  stare  mode  image  (right)  read  out  from  a  single  amplifier  and  a  corresponding 
~  10  second  chunk  of  high  frequency  data  (left).  A  few  bright  stars  and  the  corresponding 
streaks  are  indicated  by  letters.  White  arrows  in  the  left  panel  point  to  a  distinctive 
synchronous  displacement  feature  in  the  data,  visible  clearly  in  three  of  the  streaks,  in 
order  to  focus  the  reader’s  attention  to  the  non-parallelism  of  simultaneous  features  in  our 
data,  which  is  due  to  offset  in  the  original  y  position. 

While  we  achieved  significant  noise  reduction  with  our  de-trending,  our 
SNR  is  typically  a  factor  of  two  to  three  lower  than  the  Poisson-limit.  Our  noise 
is  characterized  by  high  kurtosis,  which  is  indicative  of  non-Gaussianity.  Residual 
low  frequency  fluctuations  (about  100  points  period)  are  still  noticeable  in  many  of 
our  time  series  (see  Figure  2.6).  Possible  improvements  are  discussed  in  Section  6.1. 
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Figure  2.9:  Efficiency  plotted  against  the  number  of  unidentified  candidates  (mostly  false 
positives).  Each  line  represents  a  different  window  size  re,  and  each  point  represents  the 
value  of  the  efficiency  at  threshold  9  =  0.08,  0.10,  0.15,  0.20,  0.30,  the  number  of  false 
positives  monotonically  grows  with  decreasing  9  (larger  values  of  9  on  the  left).  All  lines 
(all  w  values)  show  a  plateau  at  different  thresholds. 

2.5  Search  for  events  and  efficiency 

2.5.1  Detection  algorithm 

The  signature  of  an  occultation,  sampled  at  any  rate  >  20  Hz,  is  very 
distinctive:  it  shows  several  fluctuations  prior  to  the  Airy  ring  peak,  then  a  deep 
trough  and  possibly  a  Poisson  spot  feature,  followed  by  a  second  Airy  ring  rise  and 
more  fluctuations  (see  Figure  1.5).  The  prominence  of  these  features  depends  upon 
the  magnitude  and  spectral  type  of  the  background  star,  which  together  determine 
the  angular  size,  as  well  as  the  size  and  the  sphericity  of  the  occulter,  distance  to 
the  occulter,  and  impact  parameter  (Nihei  et  ah,  2007). 

One  possible  approach  to  detecting  occultations  in  our  lightcurves  is  to  take 
advantage  of  this  peculiar  shape,  for  example  using  correlation  of  templates,  as  in 
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Bickerton  et  al.  (2008).  Given  the  size  of  our  dataset,  however,  we  chose  to  utilize  a 
search  algorithm  general  enough  to  capture  any  fluctuation  of  some  signihcance,  but 
which  requires  less  computational  power.  We  scan  our  time  series  for  any  fluctuation 
lasting  longer  than  a  duration  w,  and  on  average  greater  than  a  threshold  9  from 
the  local  mean,  which  is  calculated  over  a  window  W  of  300  data  points  surrounding 
w.  Windows  w  of  11,  21,  31,  41  and  61  points  were  considered,  in  combination  with 
thresholds  of  0.10,  0.15,  0.20,  0.30.  We  dehne  the  local  intensity  I\{i)  as  the  ratio  of 
the  flux  in  the  local  window  w  and  in  the  surrounding  window  W .  If  the  flux  in  w 
is  suppressed  by  more  than  our  threshold  9  from  the  local  mean  (mean  over  W), 

i+w/2 

j=i—wl2 
i-\-W/2 

E  A/’-'-' 

j=j+lV/2 

then  w  is  considered  as  a  candidate.  This  is  similar  to  the  Equivalent  Width  algo¬ 
rithm,  which  is  used  in  spectral  analysis,  and  for  rare  event  searches  by  Roques  et  al. 
(2006)  and  Wang  et  al.  (2009).  Overlapping  candidates  are  then  removed  and  the 
center  of  the  window  w  that  displayed  the  largest  deviation  is  selected  as  a  single 
candidate  event.  Note  that  this  algorithm  would  in  most  cases  trigger  two  separate 
events  for  the  two  halves  of  an  occultation  on  opposite  sides  of  the  Poisson  spot 
(Figure  1.5).  These  cases  are  later  automatically  recognized  and  accounted  for  as  a 
single  event.  Different  choices  of  w  and  9  will  produce  different  detection  efficiency 
and  false  positive  rate.  We  select  an  optimized  subset  of  combinations  of  w  and  9  to 
be  used  for  our  event  detection.  This  optimization  is  described  in  the  next  section. 


<1-9, 


(2.10) 


/i(*) 
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2.5.2  Efficiency 


We  test  the  efficiency  of  our  search  by  implanting  simulated  occultations 
in  our  raw  lightcurves.  By  using  our  true  dataset  instead  of  generating  synthetic 
data  we  do  not  introduce  any  assumptions  about  the  nature  of  our  time  series.  We 
run  the  implanted  lightcurves  through  the  same  pipeline  as  the  original  lightcurves: 
de-trending  them  and  searching  for  significant  deviations  from  the  mean  flux.  In 
order  to  achieve  better  sampling  of  our  efficiency  the  entire  dataset  was  implanted 
with  one  occultation  per  lightcurve  at  each  KBO  size  we  tested: 

D  =  0.5,  0.6,  0.7,  0.8,  0.9, 1.0, 1.3,  2.0, 3.0  km,  (2.11) 


and  the  efficiency  was  assessed  for  each  size  separately.  The  finite  PSF  width  of  the 
star  induces  correlations  among  consecutive  time  stamps.  Given  the  typical  PSF  size 
in  our  data  (see  Figure  2.5)  measurements  are  considered  independent  if  separated 
by  more  than  about  seven  pixels.  Therefore,  to  modulate  the  original  time  series  by 
the  occultation  signal  we  multiply  the  star  flux  by  a  synthetic  occultation  lightcurve 
sampled  at  30  Hz.^ 

For  the  purpose  of  our  efficiency  simulations  we  assume  all  objects  are  at 
40  AU,  since  we  expect  our  occultations  to  be  within  the  Kuiper  Belt.  There  is 
little  difference  in  the  diffraction  feature  between  35  and  50  AU.  The  differences 
in  spectral  power  between  the  star  types  do  not  impact  the  occultation  features  as 

^  Since  the  occultation  typically  suppresses  the  flux,  multiplying  by  the  occultation  signal  re¬ 
duces  the  noise  by  a  factor  proportional  to  the  occultation  flux  decrease,  causing  us  to  overly 
suppress  the  Poisson  noise  by  a  factor  of  the  square  root  of  the  modulation.  Furthermore,  sources 
of  noise  that  are  not  proportional  to  the  photon  counts  (such  as  sky  background  and  read  noise) 
should  remain  constant  during  an  occultation  event,  but  this  noise  is  reduced  by  a  factor  of  the 
flux  modulation  when  the  event  is  added  to  the  lightcurve  in  this  way.  However,  since  we  expect 
to  have  a  very  high  recovery  efficiency  for  any  occultation  which  generates  effects  >  20%,  where 
the  underestimation  would  become  significant,  we  do  not  expect  this  effect  to  impact  our  efficiency 
estimation. 
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observed  by  our  system,  so  we  simulate  all  of  our  occultations  assuming  an  FOV  type 
star.  The  angular  size  of  the  star  affects  the  shape  of  the  occultation  by  smoothing 
the  diffraction  features.  It  is  therefore  important  to  properly  sample  the  angular 
size  space.  We  hnd  that,  given  the  objects  in  our  helds,  imposing  a  flat  prior  to 
the  magnitude  distribution  between  =  8  and  V  =  11  adequately  samples  our 
angular  size  range.  The  flat  prior  slightly  overestimates  the  average  cross  section 
H  of  the  events,  but  this  effect  is  more  than  compensated  by  the  loss  in  efficiency 
due  to  the  fact  that,  for  stars  with  larger  angular  sizes,  the  occultation  signal  is 
smoothed  out  as  the  diffraction  pattern  is  averaged  over  the  surface  of  the  star, 
making  the  event  harder  to  detect  (Nihei  et  ah,  2007).  Overall  our  estimate  of  our 
detection  rate  is  conservative. 

To  characterize  our  efficiency  we  implant  occultations  at  random  impact 
parameters  b  G  [0,H{D)/2].  However,  we  hrst  want  to  choose  the  most  appropriate 
window  size  and  threshold  combinations,  and  for  that  we  implant  occultations  by 
D  =  1  km  KBOs  in  the  reduced  impact  parameter  space  b  G  [0,  0.3  ■  H{D)]. 
This  set  of  modihed  lightcurves  is  used  to  optimize  our  parameters  to  maximize  our 
efficiency  and  minimize  the  number  of  false  positives  simultaneously.  Although  our 
generic  detection  approach  can  reach  high  efficiency  (nearly  100%  for  1  km  KBOs  at 
zero  impact  parameter),  it  also  produces  a  large  number  of  candidates,  most  of  which 
are  expected  to  be  false  positives.  The  combination  of  w  and  6  values  generated 
efficiencies  ranging  between  94%  (at  w=ll  and  9  =  0.1)  and  0  (at  w=Ql  and  9  =  0.3) 
and  the  number  of  candidates  ranged  between  0  and  over  1000.  Figure  2.9  shows  the 
behavior  of  the  efficiency  as  a  function  of  number  of  candidates.  Different  window 
sizes  are  represented  by  different  lines  and  the  different  thresholds  are  marked  by 
the  points  along  each  line.  Typically,  after  a  rapid  increase  in  efficiency  with  the 
decreasing  threshold,  the  efficiency  stabilizes,  while  the  nnmber  of  candidates  keeps 
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growing:  we  want  to  choose  our  parameters  near  this  point,  where  the  efficiency  is 
highest  and  any  less  stringent  choice  would  only  increase  the  number  of  our  false 
positives.  We  select  combinations  of  w  and  9  that  yield  both  an  efficiency  >  50%  and 
a  ratio  of  efficiency  to  candidates  <  0.5.  The  following  are  the  accepted  windows- 
threshold  combinations:  {w,  6)  =  (21,  0.15),  (31,  0.20),  and  (11,  0.25).  Events 
found  in  any  run  with  these  selection  parameters  were  considered  as  candidates. 
We  reached  an  overall  efficiency  of  82%  a.t  D  =  1  km  for  lightcurves  implanted  with 
synthetic  occultations  at  varying  impact  parameters  between  0  and  H /2. 

The  efficiency  of  our  search  is  summarized  in  the  top  panel  of  Figure  2.10, 
as  a  function  of  KBO  size.  We  also  plot  the  corresponding  effective  solid  angle 
r2e(-D),  dehned  as: 

fie(£')  =  E  ^  T.  6(B,  9.),  (2.12) 

* 

where  H{D,  6^)  is  the  cross  section  of  the  event,  which  depends  on  both  the  diameter 
of  the  KBO  and  the  star  angular  size  as  indicated  by  6**;  Urei  is  the  relative  velocity 
of  the  KBO,  which  depends  on  the  elongation  angle  which  is  close  to  opposition  for 
all  of  our  observations  (relative  velocities  ranging  between  24  and  25  km  s“^);  A  is 
the  distance  to  the  occulter  (assumed  to  be  A  =  40  AU),  T*  the  exposure  for  the 
star  target  (duration  of  the  lightcurve),  and  e{D,  6**)  the  recovery  efficiency  for  that 
diameter:  e{D,  6**)  =  1  if  the  implanted  event  was  recovered,  0  otherwise.  The 
sum  is  carried  out  over  all  of  our  lightcurves  with  SNR  >  25.  Qq{D)  represents  the 
equivalent  sky  coverage  of  our  survey  for  targets  at  diameter  D,  accounting  for  a 
partial  efficiency.  The  center  panel  of  Figure  2.10  shows  the  effective  solid  angle  as 
a  function  of  diameter.  The  bottom  panel  shows  the  effective  solid  angle  multiplied 
by  bracketing  slopes  for  the  size  distribution:  D~^  and  and  it  indicates  the 
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Figure  2.10:  Efficiency  of  our  survey  as  a  function  of  KBO  diameter  (top).  The  central 
panel  shows  the  effective  solid  angle  of  our  survey  and  the  bottom  panel  the  effective  solid 
angle  multiplied  D~^  and  D~*. 

survey  expects  to  see  the  largest  number  of  detections  near  D  =  700  m. 


2.5.3  Rejection  of  false  positives 

At  this  stage  we  have  more  than  a  thousand  candidates.  However,  most  of 
the  false  positives  can  be  removed  in  an  automated  fashion:  we  reject  fluctuations 
that  appear  simultaneously  in  more  than  one  lightcurve;  those  are  most  likely  due 
to  local  weather  or  atmospheric  patterns  that  were  not  corrected  in  the  de-trending 
phase  because  they  only  affected  a  subset  of  lightcurves,  and  can  be  ruled  out  as 
serendipitous  occultations.  We  also  reject  any  fluctuation  that  does  not  have  the 
right  combination  of  depth  and  width.  We  empirically  investigate  the  relationship 
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between  the  depth  and  the  width  of  an  occultation  by  a  KBO,  as  it  is  seen  by  our 
system,  taking  advantage  of  our  simulations.  To  dehne  the  depth  and  width  of  the 
events  we  fit  synthetic  occultation  lightcurves  with  inverted  top-hat  functions  with 
parameters  6  (depth)  and  7  (width).  Figure  2.11  shows  the  best  £t  values  6  and 
the  7  for  occultations  simulated  in  the  diameter  range  V  =  0.1  km  to  V  =  3.0  km, 
impact  parameters  6  =  0  to  hf/2  and  magnitude  range  8  to  11  for  FOV  stars  (the 
same  set  that  we  used  for  our  implantation  with  additional  occultations  from  objects 
D  <  0.5  km).  The  shaded  region  represents  the  area  of  this  phase  space  where  at 
least  one  occultation  was  best  fit  by  parameter  values  S  and  7  (and  the  intensity  of 
the  shade  reflects  the  frequency  of  S—j  best  fits).  We  can  automatically  reject  events 
outside  the  dashed  polygon  as  incompatible  with  D  <  3.0  km  KBO  occultations.® 
We  are  not  sensitive  to  events  shallower  than  a  10%  flux  drop. 

At  this  point  the  absolute  number  of  candidates  is  small  (25).  The  remain¬ 
ing  candidates  are  inspected  visually  (using  DS9,  Joye  &  Mandel  2003),  and  the 
lightcurves  are  extracted  with  a  different  photometric  method  (based  on  IRAF).  All 
remaining  candidates  prove  to  be  artifacts,  mostly  due  to  photometry.  No  candi¬ 
dates  are  left  after  this  elimination  process. 

2.6  Upper  limit  to  the  size  distribution  of  KBOs 
and  scientific  interpretation 

We  now  compare  our  limit  to  the  size  distribution  of  KBOs  to  that  of 
Bickerton  et  al.  (2008).  Bickerton  et  al.  (2008)  derived  an  upper  limit  to  the  surface 
density  of  KBOs  of  diameter  D  >  1  km.  They  considered  the  data  obtained  by  their 

®Note  that  the  duration  regime  over  which  we  recover  events  extend  as  far  out  as  our  largest 
window:  IF  =  61  points  or  300  ms. 
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Figure  2.11:  Phase  space  plot  showing  the  regions  of  the  flux  decrease-duration  space 
occupied  by  occultations  by  KBOs  of  diameter  0.1  to  3  km,  as  observed  through  the 
MMT/Megacam  system  bandpass.  We  simulated  occultations  from  KBOs  in  the  size 
regime  0. 1-3.0  km,  and  we  fit  the  occultations  with  an  inverted  top-hat  function  with 
parameters  5  and  7.  The  intensity  of  the  gray  scale  reflects  the  number  of  simulated 
occultations  with  best  fit  value  <5  and  7:  white  areas  are  void  of  occultations. 

own  survey  together  with  the  data  published  by  Roques  et  al.  (2006)  and  Chang 
et  al.  (2007),  assuming  100%  efficiency  for  each  survey  at  1  km,  and  obtaining  a 
total  effective  coverage  hie  =  5.4  x  10“^°  deg^.  The  cross  section  H  used  to  calculate 
hie  is  set  to  validate  the  100%  efficiency  assumption  on  a  survey  by  survey  basis. 
Our  survey  adds  7.0  x  10“®  deg^  to  the  collective  Oe,  allowing  us  to  derive  a  limit 
over  an  order  of  magnitude  stronger  than  the  limit  set  by  Bickerton  et  al.  (2008). 
Thus  we  set  a  comprehensive  95%  confidence  level  upper  limit  on  the  surface  density 
of  D  >  1  km  KBOs  at  40  AU  of  SAr(i7  >  1  km)  ~  2.0  x  10®  deg“^. 

We  can  also  derive  a  new  upper  limit  for  objects  as  small  as  700  m,  where 
our  efficiency  is  e  (  D  =  700  m)  ~  10%.  We  can  set  a  95%  conhdence  upper  limit 
of  >  0.7  km)  ~  4.8  x  10®  deg“^.  These  limits  are  shown  in  Figure  2.12, 

along  with  the  TAOS  mo  del- dependent  upper  limit  and  the  limit  set  by  the  RXTE 
X-ray  survey. 
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Figure  2.12:  Upper  limits  to  the  surface  density  of  KBOs.  Left  panel.  The  dashed  line 
is  the  best  fit  to  the  Bernstein  et  al.  (2004)  survey,  extrapolated  to  D  =  0.01km.  Three 
data-points  reported  by  Bernstein  et  al.  (2004)  are  plotted  (HST,  the  faintest  data  from 
direct  observations).  The  straight  line  is  the  TAOS  upper  limit  to  the  slope  of  the  small 
size  end  size  distribution:  q  <  4.6  (Zhang  et  ah,  2008).  The  result  by  Bickerton  et  al. 
(2008)  is  shown  as  an  empty  circle  (BKW)  as  well  as  the  X-ray  result  (Jones  et  ah,  2008, 
RXTE).  The  upper  limits  set  by  our  survey  at  U  =  1  km  and  D  =  0.7  km  are  shown  as 
filled  circles.  The  region  relevant  to  our  limit,  enclosed  in  the  square,  is  magnified  on  the 
right  hand  panel.  Right  panel.  Upper  limits  to  the  surface  density  of  KBOs,  zooming  in 
the  0.2-6  km  region  of  the  size  spectrum  where  our  survey  can  place  limits.  Symbols  and 
labels  are  the  same  as  for  the  left  panel.  The  lower  limits  for  JFC  precursor  populations 
are  also  shown  (Levison  &  Duncan,  1997;  Morbidelli,  1997;  Volk  &  Malhotra,  2008). 

2.6.1  Comparison  with  the  results  from  the  TAOS  survey 

Our  survey  aspires  to  be  complementary  to  TAOS  in  that  it  potentially 
could  detect  objects  as  small  as  300  m.  However,  at  this  stage  of  our  work  we  are 
unable  to  push  the  detection  limit  below  the  TAOS  sensitivity  (500  m).  Note  that 
the  recovery  efficiency  for  TAOS  at  700  m  is  ctaos  ~  0.3%,  a  factor  of  four  lower 
than  our  efficiency. 

The  TAOS  upper  limit  to  the  surface  density  of  KBOs  is  presented  as  a 
model-dependent  limit,  under  the  assumption  of  a  straight  power- low  behavior  for 
the  small  end  of  the  Kuiper  belt  size  distribution;  it  is  therefore  not  trivial  to  relate 
the  two  results,  but  it  is  clear  that  the  number  of  star-hours  a  dedicated  survey  can 
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collect  compensates  for  the  loss  in  efficiency  at  the  small  size  end,  and  TAOS  is  able 
to  produce  more  stringent  limits  than  our  own.  Our  survey  would  however  capture 
the  details  of  the  diffraction  feature  with  exquisite  sampling,  while  the  information 
contained  in  the  same  occultation,  as  observed  by  TAOS,  would  be  greatly  reduced 
due  to  the  slower  sampling.  This  would  allow  us  to  set  constraints  on  the  size  and 
distance  of  the  occulter,  while  the  size-distance-impact  parameter  space  is  highly 
degenerate  in  the  TAOS  data. 

2.6.2  The  Kuiper  belt  as  reservoir  of  Jupiter  Family  Comets 

The  classical  Kuiper  belt,  the  scattered  disk  objects  and  the  Plutinos  have 
all  been  considered  in  dynamical  simulations  as  possible  reservoir  of  JFCs  (see  Volk 
&  Malhotra  2008  and  references  therein).  The  inclination  distribution  of  the  JFCs 
strongly  suggests  a  disk-like  progenitor  population,  favoring  the  Kuiper  belt  over 
the  Oort  Cloud.  Giant  planets  generate  long  term  gravitational  perturbations  that 
causes  weak  orbital  chaos,  which  explains  the  injection  of  comets  to  the  JFCs  region 
(Holman  &  Wisdom,  1993;  Duncan  et  ah,  1995;  Levison  &  Duncan,  1997).  The 
efficiency  of  this  process  depends  on  the  dynamical  characteristics  of  the  progenitor 
family. 

Simulations  of  the  injection  process  lead  to  lower  limits  on  the  number  of 
progenitors,  which  we  can  compare  with  our  upper  limit  to  the  surface  density  of 
KBOs.  Bernstein  et  ah  (2004)  discussed  constraints  on  the  progenitors  of  the  JFCs 
on  the  basis  of  their  HST/ACS  survey.  This  survey  is  however  only  sensitive  to 
objects  greater  than  ~  20  km  in  diameter,  while  the  precursors  of  the  JFCs  are 
likely  to  be  in  the  size  range  1  —  10  km.  This  is  the  typical  observed  size  of  JFCs 
(Lowry  et  ah,  2008)  and  it  is  likely  that  its  progenitor  population  would  consist  of 
objects  of  similar  size  (or  slightly  larger)  than  the  JFCs  themselves. 
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In  Figure  2.12  we  show  the  lower  limits  to  the  KBO  populations  (classical 
belt  and  Plutinos)  and  scattered  disk  derived  from  dynamical  simulations.  We  use 
the  estimate  of  Levison  &  Duncan  (1997)  for  a  population  of  cometary  precursors 
entirely  in  the  classical  Kuiper  belt,  of  Morbidelli  (1997)  for  Plutinos  progenitors, 
and  of  Volk  &  Malhotra  (2008)  for  a  progenitor  population  in  the  scattered  disk. 
As  in  Bernstein  et  al.  (2004)  we  convert  the  population  estimates  for  the  Kuiper 
belt  populations  into  a  surface  density  by  assuming  for  each  population  a  projected 
sky  area  of  10^  deg^.  Volk  &  Malhotra  (2008)  provide  information  on  the  fraction 
of  time  the  objects  in  their  simulation  spend  between  30  and  50  AU  and  within  3° 
of  the  ecliptic  plane,  and  these  fractions  are  used  to  calculate  the  minimum  surface 
density  of  scattered  disk  objects  expected  in  the  region  of  sky  typically  observed  by 
occultation  surveys. 

We  are  not  presently  able  to  exclude  any  of  these  populations  as  progenitors 
of  the  JFCs.  Future  occultation  surveys,  with  improved  sensitivity,  should  provide 
valuable  information  on  the  origin  of  JFCs. 
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Chapter  3 


The  TAOS  survey  and  the 
3.75-year  dataset 


We  have  analyzed  the  first  3.75  years  of  data  from  TAOS,  the  Taiwanese 
American  Occnltation  Snrvey.  TAOS  monitors  bright  stars  to  search  for  occnltations 
by  Kniper  Belt  Objects  (KBOs).  This  dataset  comprises  5  x  10®  star-honrs  of  mnlti- 
telescope  photometric  data  taken  at  4  or  5  Hz.  No  events  consistent  with  KBO 
occnltations  were  fonnd  in  this  dataseth 


3.1  Introduction 


The  Taiwanese  American  Occnltation  Survey  (TAOS)  has  been  operating 
since  2005  with  two,  three,  and  now  four  telescopes  simultaneously  taking  stellar 
photometry  at  5  Hz^.  The  analysis  of  the  hrst  two  years  of  TAOS  data  was  described 

^This  work,  together  with  Chapter  4,  was  accepted  for  publication  by  the  Astronomical  Journal 
under  the  title  The  TAOS  Project:  Upper  Bounds  on  the  Population  of  Small  KBOs  and  Tests 
of  Models  of  Formation  and  Evolution  of  the  Outer  Solar  System,  and  it  is  the  product  of  a 
collaboration  with  the  TAOS  team  and  with  M.  J.  Holman  (Harvard-Smithsonian  Center  for 
Astrophysics) 

^A  small  subset  of  early  data  was  collected  at  4  Hz  cadence,  comprising  about  5%  of  the  data 
analyzed  in  this  work. 
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in  Zhang  et  al.  (2008,  hereinafter  Z08)  and  it  reported  no  detections.  An  upper 
limit  was  derived  to  the  slope  of  the  small  size  end  of  the  KBO  size  spectrum.  The 
TAOS  system  is  described  in  detail  in  Lehner  et  al.  (2009b).  Using  50  cm  aperture 
robotic  telescopes  in  simultaneous  observations  and  observing  with  the  relatively 
low  cadence  compared  to  all  other  KBO  occultation  surveys  (Sectionl.4),  TAOS 
was  designed  to  address  the  km-size  region  of  the  KBO  size  spectrum.  TAOS  has 
collected  over  99%  of  the  occultation  data  in  the  literature,  and  we  will  show  that 
the  marginal  sensitivity  to  sub-km  objects  is  more  than  compensated  by  the  very 
large  exposure  of  our  star  targets.  Here  we  consider  the  first  3.75  years  of  TAOS 
data,  a  signihcantly  larger  dataset  than  the  one  explored  in  Z08.  With  these  data 
we  are  able  to  begin  constraining  Kuiper  belt  formation  and  evolution  models. 

In  Section  3.2  I  describe  the  new  dataset  and  our  analysis.  In  Section  3.3 
I  briefly  describe  our  photometry  and  detection  algorithms,  as  well  as  our  efficiency 
analysis.  I  also  discuss  our  recovery  efficiency  as  a  function  of  observational  param¬ 
eters,  our  observing  strategy,  and  I  discuss  the  most  productive  strategies  for  TAOS 
and  the  other  occultation  surveys.  In  the  next  chapter  I  will  describe  the  limits 
that  this  survey  is  able  to  place  on  the  outer  solar  system  and  on  its  formation  and 
evolution  models. 


3.2  3.75  years  of  TAOS  data 

TAOS  is  a  dedicated  survey  that  observes  at  a  cadence  of  5  Hz.  The 
primary  scientihc  goal  of  the  survey  is  to  estimate  or  set  constraints  on  the  number 
of  KBOs  in  the  region  of  the  size  spectrum  that  is  currently  too  small  to  be  observed 
directly:  D  <  10  km. 
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Here  we  present  an  expanded  analysis  of  three-telescope  TAOS  data^. 
These  data  consist  of  photometric  measurements  of  target  star  helds  collected  syn¬ 
chronously  with  all  three  telescopes.  The  dataset  analyzed  here  was  collected  be¬ 
tween  January  2005  and  August  2008.  In  a  previous  analysis  of  a  subset  of  these 
data,  Z08  reported  an  upper  limit  to  the  size  distribution  of  KBOs  under  the  as¬ 
sumption  of  a  single  power  law  for  small  KBOs.  If  one  models  the  size  distribution 
for  objects  smaller  than  D  =  28  km,  the  smallest  direct  observation  (Bernstein 
et  ah,  2004),  as  a  single  power  law  dn/dD  oc  D~'^,  where  n  is  the  surface  density 
of  objects,  the  slope  of  the  distribution  is  limited  to  g  <  4.6. 

Throughout  the  remainder  of  this  paper,  a  data  run  refers  to  a  set  of  data 
collected  in  a  uninterrupted  observation  of  any  field.  For  a  single  star  in  the  held  a  set 
of  three  lightcurves  belonging  to  one  data  run  will  be  referred  to  as  a  lightcurve  set, 
and  each  three-telescope  measurement,  at  a  single  time  point,  will  be  referred  to  as  a 
triplet.  A  star-hour  refers  to  an  hour  of  high-cadence,  multi-telescope  observations 
on  a  single  target  star. 

The  data  set  described  in  this  paper  amounts  to  amount  to  5.0  x  10^  star- 
hours,  while  the  data  set  used  in  Z08  comprises  of  1.5  x  10^  star-hours.  The  details  of 
this  dataset,  and  of  the  dataset  published  in  Z08,  are  summarized  in  Table  3.1.  Over 
90%  of  our  data  is  collected  within  10°  of  the  ecliptic  plane  in  order  to  maximize 
the  rate  of  occultations. 

TAOS  uses  the  zipper-mode  technique  to  read  out  the  CCD  cameras  at 
high  frequency.  This  method,  described  in  Lehner  et  ah  (2009b),  enables  high  speed 
observations  across  the  30°  held  of  view  of  the  TAOS  telescopes,  but  it  artihcially 
increases  the  crowding  of  the  held  and  the  background.  In  zipper-mode  readout 

^The  fourth  telescope,  TAOS  C,  became  operational  in  August  2008.  The  results  presented  in 
this  paper  are  based  on  analysis  of  all  of  the  three-telescope  data  collected  to  this  point. 
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each  star  in  the  held  is  represented  in  a  subsection  of  the  output  image  -  which  we 
call  rowblock  and  which  comprises  76  rows  for  our  5  Hz  data  -  so  that  the  held  of 
view  is  entirely  imaged  in  each  rowblock.  Note  that  the  images  of  diherent  stars  in  a 
rowblock,  however,  do  not  necessarily  belong  to  the  same  time-stamp.  The  zipper¬ 
mode  readout  boosts  the  sky  background  by  a  factor  of  27  at  a  5  Hz  readout  rate. 
This  limits  the  sensitivity  of  TAOS  to  stars  as  faint  as  Mtaos  =  13.5,  for  which  a 
signal-to-noise  ratio  (SNR)  of  ~  7  can  be  achieved  in  a  dark  night.  The  magnitude 
and  SNR  distributions  for  the  target  stars  in  our  survey  are  shown  in  Figure  3.1. 
On  the  left  panel,  the  x— axis  is  the  TAOS  instrumental  magnitude  Mtaos,  which 
is  dehned  by  a  regression  on  the  USNO-B  magnitudes  to  be  similar  to  Rusno-  The 
correlation  between  instrumental  magnitude  and  SNR  is  shown  in  Figure  3.2.  The 
scatter  in  the  relationship  between  SNR  and  Mtaos  is  due  to  both  changes  in  the 
sky  background  and  in  the  weather  conditions,  and  to  different  degrees  of  crowding 
in  the  helds.  In  Figure  3.3  we  show  the  number  of  star-hours  at  different  angles  from 
opposition.  The  top  scale  indicates  the  velocity  of  a  KBO  at  the  center  of  a  held 
at  this  elongation.  We  cover  a  large  range  of  opposition  angles;  our  held  selection 
algorithm  favors  ecliptic  helds  near  zenith.  Most  angles  are  positive  because  the 
weather  at  the  site  tend  to  improve  after  midnight.  The  ehects  of  the  angle  from 
opposition  on  our  efficiency  and  event  rate,  as  well  as  the  efficiency  as  a  function  of 
magnitude  and  crowding  are  discussed  further  in  Section  3.3.5. 

3.3  Analysis 

In  this  section  we  describe  the  steps  of  the  data  analysis  procedure,  from 
photometry  to  efficiency  simulations.  In  Section  3.3.1  we  briehy  describe  the  pho¬ 
tometry  package  we  developed  to  obtain  lightcurves  from  our  CCD  images;  in  Sec- 
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Table  3.1.  TAOS  dataset  parameters  (3-telescope  data) 


Z08 

this  work 

Start  Date 

2005  February  7 

2005  February  7 

End  Date 

2006  December  31 

2008  August  2 

Light-curve  Triplets 

110,554 

366,083 

Exposure  (star-hours) 

152,787 

500,339 

Measurements^ 

7.8  X  10^ 

2.7  X  10^° 

^'Individual  telescope  measurements. 


Mfjog  star  SNR 

Figure  3.1:  Distribution  of  magnitudes  for  the  TAOS  target  stars,  (bin  size  0.16,  left). 
SNR  for  the  TAOS  target  stars,  averaged  over  the  duration  of  a  run  and  over  the  three 
telescopes  (bin  size  0.73,  right).  A  few  targets  at  greater  SNR  and  brighter  magnitude, 
amounting  to  <  5%  of  the  data,  are  not  shown. 
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Figure  3.2:  SNR  versus  TAOS  instrumental  magnitude  Mtaos!  only  a  random  sample  of 
1%  of  all  stars  is  shown  for  clarity. 
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Figure  3.3:  Distribution  of  angles  from  opposition  for  the  TAOS  targets.  The  top  axis 
shows  the  relative  velocity  of  a  KBO  at  43  AU,  given  the  position  of  the  field.  The  bin 
size  is  10.5°. 
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tion  3.3.2  we  describe  our  detection  algorithm  and  the  rejection  of  false  positives. 
We  then  describe  the  efficiency  tests:  in  Section  3.3.3  we  show  how  we  simulate 
occulation  events  to  implant  in  our  lightcurves,  and  in  Section  3.3.4  we  describe 
in  detail  how  we  identify  the  angular  size  of  our  target  stars  to  simulate  occulta- 
tions  correctly.  In  Section  3.3.5  we  describe  how  we  implant  the  occultations  in 
our  lightcurves  and  test  the  efficiency  parameters,  and  hnally  in  Section  3.3.6  we 
analyze  the  behavior  of  our  efficiency  as  a  function  of  various  parameters  relative  to 
the  occultation  and  the  observing  strategy. 

3.3.1  Photometry 

In  order  to  detect  occulations  we  need  to  identify  brief  flux  changes  in  a  star 
simultaneously  observed  by  all  telescopes.  A  custom  aperture  photometry  package 
has  been  developed  in  order  to  analyze  the  zipper-mode  TAOS  images  (Zhang  et  ah, 
2009).  TAOS  photometry  is  conducted  in  three  phases:  star  identihcation  and 
magnitude  assignment,  determination  of  the  aperture  size,  and  photometry. 

In  order  to  identify  the  star  targets  we  collect  conventional,  or  stare-mode, 
images  of  our  target  held  at  the  beginning  of  each  observing  run.  These  images  are 
analyzed  using  SExtractor  (Bertin  &  Arnouts,  1996).  The  sources  thus  extracted 
are  identihed  with  USNO-B  catalog  sources  (Monet  et  ah,  2003)  using  the  WCS 
tools  (Mink,  2006).  A  magnitude  regression  between  USNO-B  R2  and  the  instru¬ 
mental  magnitude  is  computed  for  all  stars  with  a  USNO-B  counterpart  that  are 
brighter  than  Mtaos  =  13.5:  those  are  hereinafter  identihed  as  targets. 

Next,  the  best  aperture  size  for  each  target  star  is  empirically  determined. 
Because  of  the  problems  caused  by  the  artihcial  crowding  of  the  zipper-mode  images 
a  diherent  aperture  size  for  each  star  is  chosen  so  as  to  minimize  contamination.  For 
each  star  the  SNR  of  the  star  is  evaluated  as  a  function  of  aperture  size  for  the  hrst 
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~  1000  time-stamps  and  the  size  that  generates  the  best  SNR  is  chosen  for  that 
star.  In  this  fashion,  wherever  near-by  star  images  would  contaminate  the  target, 
the  aperture  is  naturally  chosen  conservatively.  Square  apertures  are  employed  to 
increase  the  computational  efficiency. 

The  initial  position  of  the  aperture  for  every  target  is  estimated  from  the 
stare-mode  image.  Corrections  are  made  to  the  aperture  centers  from  time-stamp 
to  time-stamp  to  account  for  small  translational  or  rotational  displacements  of  the 
held.  The  sky  is  estimated  as  the  mode  of  the  hux  in  each  column  in  each  rowblock, 
and  then  it  is  subtracted  from  the  star  hux.  Note  that  the  zipper-mode  readout 
induces  excess  counts  in  the  columns  where  bright  stars  are  located  since  the  charges 
are  transferred  with  the  shutter  open  between  each  exposure.  Estimating  the  back¬ 
ground  column  by  column  allows  us  to  remove  the  sky  background  as  well  as  the 
excess  hux  in  these  streaks.  Time-stamp  by  time-stamp  the  hux  for  each  targets  is 
thus  collected,  composing  a  photometric  time  series  for  each  target  star. 

3.3.2  Detection  and  false  positive  rejection 

The  Fresnel  scale  is  dehned  as  F  =  (AA/2)^  where  A  is  the  wavelength  of 
observation,  and  A  the  distance  to  the  occulter  (Roques  et  ah,  1987;  Born  &  Wolf, 
1980).  For  optical  observations  at  the  distance  of  the  Kuiper  belt  (about  43  AU)  the 
Fresnel  scale  is  F  ~  1.4  km.  Dihraction  ehects  are  therefore  important  in  the  events 
that  we  are  seeking.  Occultations  are  manifested  in  the  lightcurve  of  an  observed  star 
as  an  alternation  of  bright  and  dark  features,  typically  with  an  overall  suppression 
of  the  hux.  Theoretical  occultation  lightcurves  are  shown  in  Figures  3.4  &  3.5.  The 
signature  of  an  occultation  by  a  KBO  of  sub-kilometer  size  has  a  typical  duration 
of  about  0.2  second  at  opposition,  and  about  a  second  near  quadrature.  A  typical 
KBO  occultation  is  then  expected  to  result  in  the  suppression  of  the  hux  for  one 
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or  a  few  consecutive  points  in  a  TAOS  lightcurve.  In  order  to  ascertain  the  extra¬ 
terrestrial  origin  of  a  dip  in  a  lightcurve  TAOS  observes  simultaneously  with  multiple 
telescopes.  This  allows  us  to  rule  out,  on  the  basis  of  simple  parallax  considerations, 
atmospheric  scintillation  phenomena  which  might  mimic  an  occultation  event  and 
which  could  be  a  source  of  false  positives  in  occultation  surveys,  as  well  as  any 
non-atmospheric  phenomena  such  as  birds,  airplanes,  etc. 

The  statistical  signihcance  of  a  simultaneous  low  point  in  our  lightcurves 
can  be  assessed  rigorously,  and  the  probability  of  a  low  measurement  being  drawn  out 
of  pure  noise  decreases  with  the  number  of  telescopes  observing  the  target,  provided 
that  the  measurements  for  the  telescopes  are  independent.  The  lightcurves  are  high- 
pass  hltered  to  remove  trends  due  to  weather  patterns  and  changes  in  atmospheric 
transparency.  High-pass  hltering  the  lightcurves  preserves  the  information  on  time- 
scales  relevant  to  occultation  phenomena  (one  or  a  few  points  in  a  time  series).  The 
implementation  of  the  hlter  is  described  in  Z08.  The  hlter  produces  a  time  series 
in  which  h{t),  the  measurement  taken  at  time  t,  represents  the  deviation  from  the 
local  mean  of  the  lightcurve  in  units  of  local  standard  deviation.  Any  run  where  the 
independence  of  the  measurements  after  hltering  cannot  be  rigorously  established  is 
removed  from  our  dataset  (Lehner  et  ah,  2009a). 

To  detect  events  we  rank-order  the  photometric  measurement  in  each  of 
our  light-curves,  from  the  lowest  to  the  highest  hux,  independently  for  each  telescope 
(labeled  A,  B  and  D).  The  i-th  point  in  a  lightcurve  will  be  associated  to  rank  rj 
for  telescope  T.  We  then  consider  the  rank  triplets  (r)^,  rf,  rf).  The  probability 
distribution  of  the  quantity  Zi  =  —  ln{rf  rf  rf  /N^},  with  Np  the  number  of 

points  in  the  lightcurve  set,  can  be  determined  combinatorially.  Knowing  this, 
under  the  null  hypothesis  that  there  is  no  event  in  the  triplet  i,  we  can  compute  the 
probability  for  a  random  variable  Z  arising  from  this  distribution  P{Z  >  Zi)  = 
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We  set  a  threshold  such  that  we  expect  fewer  than  0.27  events  in  our  dataset  that 
are  due  to  random  fluctuations.  For  the  dataset  discussed  in  this  paper  we  accept  as 
events  all  data  points  that  produce  a  rank  product  less  likely  than  ^  =  3.0  x  10“^^ 
to  be  drawn  from  a  random  distribution.  Note  that  events  generated  by  large 
KBOs,  or  for  observations  near  quadrature,  would  affect  more  than  one  point  in 
the  lightcurve  (Figure  3.5),  and  our  rank-based  search  algorithm  is  most  efficient 
when  the  dip  in  the  lightcurve  is  isolated.  Therefore,  in  addition  to  searching  for 
single-point  events,  we  also  bin  our  lightcurves  by  2,  re-rank  them  and  repeat  the 
statistical  tests  described  above.  The  probability  of  each  data  point  is  assessed  for 
both  unbinned  and  binned  lightcurves.  Each  lightcurve  is  binned  twice,  with  two 
different  starting  points.  This  increases  the  detectability  of  occultations  by  large 
KBOs  and  by  KBOs  transiting  with  low  relative  velocity.  For  a  detailed  discussion 
of  our  statistical  analysis  see  Lehner  et  ah  (2009a). 

For  a  set  of  lightcurves  of  a  given  star,  the  ranks  in  the  three  telescopes 
should  not  be  correlated  for  the  statistical  analysis  described  above  to  be  valid.  A 
three-dimensional  scatter-plot  of  the  rank  triplets  should  thus  uniformly  populate 
a  cube.  We  have  devised  two  statistical  tests  to  assess  the  quality  of  a  data  run 
which  are  based  on  this  requirement.  First,  for  each  lightcurve  set  we  divide  the 
scatter-plot  into  a  three-dimensional  grid  and  count  how  many  rank  triplets  are  in 
each  element.  We  expect  the  number  in  each  element  to  be  equal,  and  we  thus 
compute  a  Pearsons  statistic  to  measure  how  far  the  distribution  deviates  from 
the  expected  uniform  distribution.  If  any  data  runs  show  a  signihcant  deviation 
from  the  predicted  distribution,  the  entire  data  run  is  rejected.  Second,  we  are 
particularly  concerned  with  any  correlation  which  would  produce  a  bias  toward 
simultaneously  low  ranks.  We  thus  look  at  the  data  in  the  grid  element  with  the 
lowest  ranks  in  all  telescopes.  The  number  of  rank  triplets  in  this  grid  element 
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follows  the  hypergeometric  distribution  (Rice,  2006,  pg.  39).  Once  again,  we  test  to 
see  if  the  number  of  triplets  in  this  grid  element  follows  the  predicted  distribution  for 
all  of  the  lightcurves,  and  once  again  we  reject  any  data  run  which  shows  a  signihcant 
deviation.  These  tests  will  be  described  in  detail  in  Lehner  et  ah  (2009a).  After 
removing  the  data  runs  that  do  not  fulhll  our  independence  requirement  we  have 
8  X  10^  star-hours  of  data. 

For  the  next  step  we  relax  considerably  the  ultimate  signihcance  require¬ 
ment  described  above,  and  select  as  provisional  candidates  those  triplets  that  have 
^  <  1.0  X  10“®.  Note  that  the  signihcance  ^  refers  to  the  probability  that  the  point 
would  be  drawn  from  a  random  distribution,  therefore  the  lower  the  value  of  ^  the 
higher  the  statistical  signihcance  of  the  event.  Nearly  150,000  provisional  candidates 
are  found.  We  use  all  these  measurements  to  identify  and  remove  spurious  regions  of 
the  lightcurves,  and  hence  identify  and  reject  false  positive  events  which  arise  from 
sources  other  than  random  chance.  The  constraints  described  below  allows  us  to 
recognize  regions  of  lightcurves  with  atypical  noise  and  contamination  by  transiting 
objects  (satellites  or  meteors,  which  turn  out  to  be  the  major  source  of  false  alarms), 
and  to  identify  high-frequency  huctuations  in  the  raw  data  that  are  not  removed  by 
the  high-pass  hlter.  These  are  the  steps  of  our  false  positive  rejection  process: 

•  Contiguity:  Contiguous  candidates  within  a  lightcurve  and  candidates  that 
are  within  three  time-stamps  of  each  other  are  removed.  Only  the  one  rank 
triplet  that  has  the  highest  signihcance  in  a  series  of  contiguous  or  proximate 
points  is  considered  as  a  candidate.  This  removes  double-counted  events: 
events  caused  by  large  KBOs  or  KBOs  moving  at  low  relative  velocity  would 
affect  more  than  one  contiguous  point.  Furthermore  this  removes  events  that 
are  double-counted  because  they  appear  signihcant  in  both  the  binned  and 
unbinned  lightcurves.  This  eliminates  about  40%  of  the  candidates. 
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•  Simultaneity:  Candidates  that  appear  in  the  lightcurves  of  more  than  one  star 
simultaneously  at  the  same  time-stamp  or  within  three  rowblocks  are  consid¬ 
ered  to  be  false  positives.  We  expect  simultaneous  count  drops  in  time-domain 
to  be  primarily  due  to  inaccurate  aperture  positioning  in  the  photometry.  In 
the  rowblock  domain  simultaneous  count  drops  might  be  due  to  inaccurate 
background  determination  or  non-occultation  events  altering  the  baseline  of 
the  lightcurve  at,  or  around,  the  candidate  event,  or  by  fast  moving  cirrus 
clouds  or  other  phenomena  which  induce  high  frequency  fluctuations  in  multi¬ 
ple  lightcurves  which  are  not  removed  by  the  high-pass  hlter.  This  cut  removes 
about  60%  of  the  remaining  candidates. 

•  Number  of  telescopes:  At  this  point  we  require  all  of  our  remaining  candidates 
to  have  been  observed  by  all  three  telescopes.  Although  we  are  only  considering 
three-telescope  runs  in  our  analysis,  for  some  targets  the  lightcurve  might  not 
be  extracted  in  the  photometry  phase  for  all  telescopes.  Small  differences 
in  the  held  of  view  and  in  the  held  distortion  might  make  one  star  target 
not  visible  to  all  telescopes  if  it  is  at  the  edge  of  the  held  or  if  the  crowding 
induced  by  the  zipper-mode  readout  caused  overlap  of  the  target  with  other 
stars  (Lehner  et  ah,  2009b).  About  35%  of  the  remaining  candidates  are  thus 
removed.  At  this  point  there  are  still  over  20,000  candidates  left.  Note  that 
we  do  not  count  the  star  hours  discarded  by  this  cut  in  our  total  exposure  of 
5  X  10^  star-hours. 

•  Significance  threshold:  We  hnally  constrain  ^  such  as  to  expect  fewer  than 
0.27  false  positives  in  our  dataset.  This  constraint  depends  on  the  size  of  the 
dataset:  for  the  9  x  10®  triplets  remaining  ^  <  3.0  x  10“^^  allows  <  0.27  false 
positives  due  to  random  noise.  Only  228  candidates  remain. 
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The  remaining  candidates  require  visual  inspection:  first  of  the  lightcurves, 
and  for  any  remaining  candidates,  of  the  images.  Most  of  the  events  are  caused  by 
the  passing  of  bright  objects,  such  as  artihcial  satellites,  meteorites  or  asteroids, 
that  generate  a  variation  in  the  background  or  baseline  of  the  lightcurve  responsible 
for  causing  artihcially  low  counts  in  the  neighborhood  of  the  object.  Many,  but  not 
all,  of  these  false  positives  are  removed  by  the  simultaneity  cut  described  above. 
In  the  presence  of  a  bright  object  overlapping  with  a  star-streak  generated  by  the 
zipper-mode  readout,  the  brightness  of  the  streak  is  overestimated,  thus  too  much 
flux  is  subtracted  from  the  rowblock  column  causing  an  artihcial  hux  drop  in  the 
star  time  series.  In  many  instances  the  foreground  object  will  also  appear  inside 
the  star  aperture  artihcially  boosting  its  brightness.  This  hux  drop  will  then  be 
associated  with  a  very  high  hux  measurement  following  or  preceding  the  event  epoch, 
a  signature  that  allows  us  to  remove  these  false  positives  by  inspecting  the  lightcurve. 
We  also  inspect  the  centroid  position  of  the  aperture.  If  the  aperture  position  has 
moved  signihcantly  at  the  time-stamp  of  the  candidate  the  candidate  is  rejected.  Of 
the  remaining  candidates,  90%  are  rejected  by  visual  inspection  of  the  time  series. 

Finally  we  inspect  the  images  of  the  remaining  23  candidates:  they  also 
were  all  associated  with  bright  moving  objects  overlapping  star  streaks.  No  candi¬ 
date  events  were  left  in  our  dataset  at  the  conclusion  of  this  process. 

3.3.3  Occultation  event  simulator 

It  is  necessary  to  assess  the  efficiency  of  our  recovery  algorithm  in  order 
to  derive  the  number  density  of  KBOs  from  the  number  of  events  in  our  survey.  In 
order  to  measure  our  recovery  efficiency  we  implant  our  data  with  synthetic  occul- 
tations.  By  implanting  into  the  actual  lightcurves  we  do  not  make  any  assumption 
regarding  the  nature  of  the  noise  in  our  data.  Note  that  our  detection  algorithm. 
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described  in  Section  3.3.1,  which  is  based  on  rank  statistics,  is  not  affected  by  the 
spectral  characteristics  of  the  noise,  as  long  as  the  noise  is  stationary  (Lehner  et  ah, 
2009a).  It  is  important  to  correctly  characterize  occultations  by  all  objects  we  might 
detect,  and  for  this  we  need  to  properly  characterize  all  of  our  target  stars  as  well 
as  the  effects  of  our  telescope  and  camera  system.  Here  we  describe  the  steps  of 
our  algorithms  to  generate  synthetic  occultation  features  to  be  implanted  in  our 
lightcurves.  Our  simulator  is  also  described  in  Nihei  et  ah  (2007). 

For  a  spherical  occulter,  and  under  the  plane  wave  approximation,  the 
intensity  of  a  point  on  a  diffraction  pattern  is  uniquely  determined  by  four  parame¬ 
ters:  the  wavelength  A,  the  distance  to  the  occulting  object  A,  the  diameter  of  the 
object  D,  and  the  perpendicular  distance,  x,  of  that  point  of  the  diffraction  pattern 
from  a  line  that  passes  through  the  star  and  the  occulter.  We  can  uniquely  describe 
the  intensity  at  any  point  of  an  occultation  pattern  (or  lightcurve)  as  a  function  of 
two  dimensionless  quantities  S  =  D/F{X,  A)  and  r]  =  x/F{X,  A),  where  F  is  the 
Fresnel  scale  (Roques  et  ah,  1987,  and  references  therein). 

The  calculation  of  the  intensity  pattern  requires  signihcant  computational 
time,  as  a  result  of  the  highly  oscillatory  nature  of  the  integrals  involved  in  the 
calculation  (Nihei  et  ah,  2007).  We  therefore  compiled  a  diffraction  table:  a  binary 
hie  that  contains  the  intensity  as  a  function  of  5  and  rj.  We  can  generate  a  theoretical 
diffraction  lightcurve  by  interpolating  in  this  table.  A  gray  scale  representation  of 
the  diffraction  table  is  shown  in  Figure  3.4,  along  with  an  example  of  a  diffraction 
pattern  for  a.  D  =  3.0  km  KBO.  Note  the  following  features  in  the  diffraction  pattern: 
the  Airy  peaks,  the  brightest  fringes  at  r/  =  ±2.4  [x  =  ±3.36  km),  and  the  Poisson 
spot,  the  central  brightening  (Roques  et  ah,  1987). 

We  also  need  to  simulate  our  telescope  and  camera  system.  The  TAOS  sys¬ 
tem  uses  wide  bandpass  hlters  with  an  approximately  constant  response  between  500 
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Figure  3.4:  Diffraction  table  (top),  the  y-axis  is  the  diameter  (5)  of  the  occulter  in  Fresnel 
units,  and  the  x-axis  (77)  is  the  distance  from  the  center  of  the  occultation  in  Fresnel  units. 
Physical  units  are  marked  on  the  right  and  top  axes  in  km,  assuming  KBOs  at  43  AU 
and  A  =  600  nm.  The  gray  scale  represents  the  intensity  of  the  flux.  At  the  bottom  an 
occultation  for  a  D  =  3  km  KBO,  or  (5  =  2.16,  as  extracted  from  the  table.  The  location 
of  the  light-curve  on  the  table  is  indicated  by  a  dashed  line 
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and  700  nm  (Lehner  et  al.,  2009b).  In  order  to  characterize  the  occultations  properly 
as  they  would  be  seen  by  the  TAOS  telescopes  we  weight  the  stellar  spectrum  s(A) 
by  the  TAOS  transmission  /(A)  (hlter  bandpass  and  camera  quantum  efficiency  as 
a  function  of  wavelength)  and  then  we  calculate  the  diffraction  intensity  point  by 
point: 

/OO 

(3-1) 

•OO 

as  described  in  Nihei  et  ah  (2007).  We  found  that  the  differences  in  the  spectra 
among  our  target  stars  has  a  negligible  impact  on  the  shape  of  the  occultation 
pattern  as  observed  by  our  system.  We  therefore  simulate  all  of  our  occultations 
assuming  an  FOV  spectral  type.  A  lightcurve  generated  with  our  simulator  for  a 
3  km  KBO  at  43  AU  occulting  an  FOV  point  source  is  shown  in  Figure  3.5,  top. 

We  now  have  a  diffraction  lightcurve  for  our  KBO  occulting  a  point  source. 
Most  of  our  targets,  however,  have  a  hnite  size  projected  on  the  plane  of  the  occulter 
in  the  sky.  We  have  to  take  such  hnite  source  effects  into  account,  since  integrating 
the  occultation  over  the  disk  of  the  star  will,  at  times  signihcantly,  smooth  the 
occultation  features,  while  broadening  the  occultation  pattern  (and  thus  the  event 
cross  section).  The  determination  of  the  angular  size  for  our  targets  is  discussed  in 
Section  3.3.4. 

Keeping  the  stellar  type  hxed,  the  angular  size  is  modulated  by  changing 
the  apparent  magnitude  of  the  star  and  we  can  use  the  point  source  lightcurve  to 
integrate  the  occultation  signature  over  the  star  disk.  For  every  star  and  KBO  pair 
for  which  we  generated  a  point  source  lightcurve  we  generate  hnite  source  lightcurves 
for  stars  in  the  magnitude  range  8  <  V  <14,  with  steps  of  0.1  mag.  A  magnitude 
V  =  11  hnite  source  lightcurve  generated  from  the  point  source  in  the  top  panel  of 
Figure  3.5,  is  shown  in  the  second  panel,  left.  Note  the  smoothing  of  the  dihraction 
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features. 


We  then  proceed  to  implant  occultation  events  in  our  data,  which  are 
then  reprocessed  in  the  same  way  we  did  to  search  for  true  events.  To  generate  an 
event,  we  hrst  determine  the  star  parameters  by  identifying  the  star  in  the  USNO-B 
catalog  (Monet  et  ah,  2003).  The  elongation  angle  of  the  observation  determines  the 
relative  velocity  of  the  KBO,  together  with  the  KBO  distance  (see  Equation  1.8) 
We  choose  the  parameters  for  the  occultation:  we  select  a  diameter  ac¬ 
cording  to  a  distribution  that  will  allow  us  to  have  sufficient  number  of  recoveries 
at  all  sizes,  so  as  to  be  able  to  sample  our  efficiency  properly  (this  will  be  discussed 
in  Section  3.3.5),  and  we  draw  a  random  impact  parameter  and  occultation  epoch 
in  the  lightcurve.  We  generate  the  occultation  lightcurve  and  implant  its  signal  into 
our  star  lightcurve: 

•  We  determine  the  angular  size  of  the  target  star  (see  Section  3.3.4),  match 
it  to  a  previously  generated  hnite  source  lightcurve  and  load  the  hnite  source 
lightcurve. 

•  We  modify  the  lightcurve  to  account  for  a  hnite  impact  parameter  b:  for  this 
we  use  the  hnite  source  lightcurve  at  impact  parameter  6  =  0  as  an  input, 
calculating  the  intensity  of  the  occultation  signal  at  the  new  distance  of  each 
point  by  interpolating  points  of  the  hnite  source  lightcurve.  A  lightcurve  for  an 
FOV,  V=ll  star  and  a  3  km  KBO  occulting  at  an  impact  parameter  b  =  2  km 
is  plotted  on  the  right  hand  side  of  the  second  panel  of  Figure  3.5. 

•  We  then  calculate  the  relative  velocity  of  the  KBO  as  a  function  of  distance 
and  angle  from  opposition.  Knowing  this  we  can  confer  time  information  to 
our  occultation  features,  which  were  thus  far  expressed  as  a  function  of  the 
distance  to  the  point  of  minimum  approach,  and  set  the  time-line  of  the  event. 
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Figure  3.5:  Steps  of  the  generation  of  a  simulated  occultation  event.  Top,  the  left  and  right 
panels  both  show  the  same  point  source  lightcurve  for  a  3  km  KBO  at  43  AU  occulting  an 
FOV  star.  Second  row:  finite  source  lightcurves  for  the  same  occultation  parameters  for  a 
V  =  11  star,  corresponding  to  an  angular  size  of  0.015  mas,  for  a  zero  impact  parameter 
(left),  and  at  an  impact  parameter  of  2  km  (right).  Row  three:  the  lightcurves  in  row  two 
are  integrated  over  intervals  of  105  ms  for  the  occultation  above  at  opposition  (left)  and 
at  50°  from  opposition  (right).  The  lightcurves  are  sampled  at  5  Hz,  with  no  time  offset 
(left)  and  with  a  time  offset  of  50  ms  (right). 


•  Finally  we  smooth  the  lightcurve  to  account  for  hnite  exposure  intervals,  and 
we  sample  the  finite  exposure  lightcurve  at  the  appropriate  sampling  rate.  In 
this  step  we  can  account  for  dead-time  in  the  sampling  interval,  which  for 
TAOS  is  47.5%  at  5  Hz.  We  also  allow  an  offset  in  time  between  the  center 
of  the  hnite  sampled  lightcurve  and  the  integration  bin.  In  the  bottom  row  of 
Figure  3.5  the  lightcurves  for  the  event  in  row  two  are  integrated  over  105  ms 
second  intervals  and  sampled  at  5  Hz,  the  typical  sampling  rate  of  TAOS,  for 
an  event  at  opposition  and  with  no  time  offset  (left)  and  an  event  at  50°  from 
opposition  with  an  offset  of  50  ms  between  the  center  of  the  sampling  interval 
and  the  center  of  the  occultation  (right). 


We  now  have  a  proper  occultation  signature,  which  can  be  implanted  in 


our  light-curves. 


3.3.4  Determination  of  the  stellar  angular  size 

The  angular  size  of  the  target  star  is  an  important  parameter  in  the  gen¬ 
eration  of  an  occultation  pattern.  Our  helds  contain  a  variety  of  stellar  types  and  a 
large  range  of  angular  sizes  (Figure  3.6a). 

Angular  sizes  have  been  related  to  the  position  of  a  star  in  the  color-color 
or  color- magnitude  diagrams  (e.g.,  van  Belle  1999,  Nordgren  et  al.  2002).  We  follow 
the  work  of  Nordgren  et  al.  (2002)  and  calculate  the  angular  size  of  our  star  targets 
using  the  2MASS  J  and  K  color  (Cutri  et  al.,  2003)  to  invert  the  set  of  equations: 

Fk  =  (3.942  ±  0.006)  -  (0.095  ±0.007)(J- AT)  (3.2) 

Fk  =  4.2207 -O.IA:- 0.5  log 0  (3.3) 

where  Fx  is  the  surface  brightness  of  a  star  in  K-band,  which  is  related  to  its  J  —  K 
color,  as  well  as  to  its  unreddened  apparent  K  magnitude  and  angular  size  6.  The 
relationship  between  the  surface  brightness  and  the  color  of  a  star  (Equation  3.2) 
is  calibrated  using  angular  sizes  measured  directly  by  long  baseline  interferome¬ 
try  (Nordgren  et  al.,  2002). 

Not  ah  of  our  target  stars,  however,  are  identihed  2MASS  objects,  while  in 
the  photometry  phase  we  have  identihed  all  of  our  targets  with  USNO-B  objects.  We 
therefore  devised  a  method  that  relies  on  USNO-B  R  and  B  magnitude  to  calculate 
the  angular  sizes  of  our  targets. 

We  hrst  derive  the  angular  size  of  a  subset  of  targets  identihed  with  2MASS 
objects  using  the  above  equations,  and  scale  it  to  obtain  the  angular  size  the  targets 
would  have  if  their  apparent  magnitude  were  R  =  12.  We  then  considered  the 
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USNO-B  B  —  R  color  for  all  of  these  targets  and  calculated  a  regression  on  these 
points.  This  generates  a  formula  that  allows  us  to  go  from  the  USNO-B  color  of 
any  of  our  targets  to  2MASS  colors  and  thus  predict  angular  sizes  according  to 
Equations  3.2  &  3.3,  for  an  apparent  magnitude  R  =  12.  To  calculate  the  true 
angular  size  we  rescale  from  i?  =  12  to  i?usNo"^-  The  angular  sizes  of  a  subset  of 
TAOS  targets,  rescaled  to  R  =  12,  is  plotted  as  derived  from  Equation  3.2  &  3.3 
versus  the  USNO-B  B  —  R  color  (Figure  3.6b).  Our  regression  on  the  data  is  plotted 
as  well  (solid  line). 

The  scatter  in  the  determination  of  the  angular  size  via  the  method  de¬ 
scribed  above  is  large,  as  can  be  seen  in  Figure  3.6b.  This  is  due  to  scatter  in  the 
USNO-B  color  (^  0.3  mag,  Monet  et  al.  2003),  to  the  (much  smaller)  scatter  in 
the  J  and  K  magnitudes  and  to  the  scatter  in  the  empirical  determination  of  the 
relationship  between  9  and  J  —  A'  in  Equation  3.2  &  3.3.  We  have  not  used  any 
interstellar  reddening  corrections,  and  the  angular  size  estimation  of  an  unknown 
reddened  star  from  the  near-IR  relationship  would  be  relatively  less  affected  com¬ 
pared  to  that  in  visual  bands.  Reddening  is  typically  small  for  our  targets  though, 
since  we  are  only  considering  objects  brighter  than  R  ~  13.5. 

The  distribution  of  angular  sizes  is  well  reproduced.  Figure  3.6a  shows  the 
distribution  of  angular  sizes  for  a  typical  TAOS  held,  calculated  via  Equations  3.2 
&  3.3,  and  Figure  3.6c  shows  the  distribution  of  angular  sizes  in  our  efficiency 
simulation  obtained  via  the  USNO-B  color.  The  distributions  do  overlap.  About 
2%  of  our  simulated  angular  sizes  fall  in  the  region  9  >  0.15  mas  and  Rtaos  > 
11,  where  there  are  no  observed  objects.  These  objects  have  poor  USNO-B  color 
determination.  Figure  3.6d  shows  the  region  of  the  6*  — Rtaos  space  where  simulated 

^We  do  not  use  our  instrumental  magnitude  for  rescaling  for  consistency  with  what  is  used  in 
the  color  determination. 
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Figure  3.6:  Angular  size  distribution  for  a  typical  TAOS  field  (field  120,  RA:  13°. 7,  Dec: 
—  10°. 7)  derived  from  the  2MASS  K  —  J  colors  (a).  The  curves  show  the  theoretical 
behavior  of  the  angular  size  for  AO,  FO,  GO,  KO  dwarf  stars  (thin  lines)  and  GO,  KO,  K5 
giants  (thick  lines).  The  size  of  the  Fresnel  scale  at  43  AU  is  shown  at  0.08  mas  (dashed 
line).  Best  fit  to  the  angular  size  distribution  (b):  x-axis  is  the  USNO-B  B  —  R  color  and 
y  is  the  angular  size  derived  using  Equations  3.2  &:  3.3,  converting  all  stars  to  apparent 
magnitude  R  =  12.  Implanted  angular  sizes,  derived  from  USNO-B  the  B  —  R  color 
for  our  simulation  (c)  and  angular  size  of  the  star  for  which  events  are  recovered  (d);  a 
random  subset  of  1%  of  our  data  is  plotted  in  the  bottom  panel. 


events  are  recovered.  There  are  few  recoveries  in  the  region  6  >  0.15  mas  and 
.Rtaos  >  11,  so  these  stars  do  not  contribute  the  the  expected  event  rate. 


3.3.5  Implantation  and  efficiency  test  parameters 

We  now  assess  our  efficiency  at  detecting  events  in  the  Kuiper  Belt.  In 
order  to  sample  properly  the  space  of  diameters  to  which  the  survey  is  sensitive  we 
generate  synthetic  occultations  by  objects  of  diameter  D  =  0.5,  0.7,  1.0,  1.3,  2.0,  3.0, 
8.0,  and  30.0  km.  For  a  30  km  diameter  KBO  the  event  falls  in  the  geometric  regime. 
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Table  3.2.  Distribution  of  synthetic  events 


diameter  (km) 

implantations 

recoveries 

30.0 

231 

75 

8.0 

385 

84 

3.0 

1078 

73 

2.0 

2003 

89 

1.3 

4393 

73 

1.0 

13255 

66 

0.7 

36222 

40 

0.5 

447764 

9 

diffraction  effects  are  therefore  no  longer  signihcant  and  our  efficiency  stabilizes. 
Because  our  sensitivity  decreases  with  decreasing  diameter  we  implant  progressively 
more  objects  at  smaller  diameters.  The  number  of  implantations  at  each  size  is 
designed  to  allow  us  to  obtain  a  good  sampling  at  all  sizes.  In  Table  3.2  we  report 
the  number  of  objects  implanted  for  each  size  in  one  of  our  efficiency  runs,  and  the 
number  of  recoveries^.  For  objects  within  the  Kuiper  belt  (about  30  to  60  AU),  the 
differences  induced  by  different  distances  are  negligible  in  the  occultation  features 
as  observed  by  TAOS.  We  therefore  set  the  distance  to  A  =  43  AU. 

Every  occultation  event  is  implanted  at  a  random  epoch  in  the  lightcurve 
set.  We  also  draw  a  random  impact  parameter  between  0  and  H/2,  where  we  set 
H,  a  measure  of  the  cross  section  of  the  event,  to  the  size  of  the  Airy  ring  and  the 
projected  size  of  the  star,  as  descibed  in  Section  1.3.1  (Equation  1.3.1). 

In  order  to  implant  these  occultations  we  modulate  the  lightcurve  by 
subtracting  (adding)  the  amount  of  flux  suppressed  (augmented)  by  the  occul¬ 
tation  at  each  data-point:  that  is,  the  implanted  lightcurve  focc  is  generated  as 
foccit)  =  fit)  —  (1  —  d(t))/,  where  /  is  the  local  mean  of  the  flux,  and  df)  is 
^Four  runs  are  conducted  to  reduce  the  error  in  the  determination  of  our  efficiency. 
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Figure  3.7:  Implanted  occultations  recovered  by  our  pipeline.  On  the  left  the  lightcurves 
refer  to  star  targets  in  a  field  near  opposition  (Field  158,  RA:  18*^  53™  44®,  Dec:  - 
17°  23^  00^')  and  on  the  right  to  a  held  at  large  angle  from  opposition  (Field  22,  RA: 
4h  52™  43®^  Dec:  —20°  51'  00"). Each  data-point  h{t)  represents  the  deviation  of  that 
time-stamp  from  the  local  mean  of  the  lightcurve  in  units  of  standard  deviation.  The 
lightcurve  for  telescope  D  is  plotted  in  each  panel,  while  the  relevant  measurement,  h{0), 
for  telescope  A  and  B  are  indicated  by  arrows  with  the  height  of  the  arrow  representing  the 
measurement  at  the  center  of  the  occultation  t  =  0.  Each  event  is  described  in  Table  3.3 


the  simulated  event  lightcurve  with  baseline  d{t)  1  away  from  the  event.  This 
approach  slightly  overestimates  the  noise  due  to  Poisson  statistics  where  the  flux  is 
suppressed,  giving  us  a  conservative  estimate  of  our  efficiency.  We  implant  exactly 
one  occultation  in  each  lightcurve  in  our  dataset. 


3.3.6  Analysis  of  the  efficiency  parameters 

We  now  process  onr  implanted  lightcurves  as  we  previously  did  to  search 
for  events,  namely  we  hlter,  rank  the  lightcnrves,  and  evaluate  the  signihcance  ^  of 
each  point  in  each  lightcnrve.  We  then  remove  the  false  positives  as  described  in 
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Table  3.3.  Parameter  of  implanted  events  in  Figure  3.7 


77  (km) 

b  (km) 

Urel  (km/s) 

SNR 

6**  (mas) 

(a) 

0.7 

0.50 

25.4 

39.8 

0.03 

(b) 

1.0 

0.50 

14.9 

29.4 

0.03 

(c) 

3.0 

1.19 

25.4 

12.8 

0.03 

(d) 

8.0 

2.7 

25.4 

9.4 

0.006 

(e) 

0.7 

0.22 

8.2 

10.2 

0.03 

(f) 

1.0 

0.94 

7.7 

10.8 

0.03 

(g) 

3.0 

1.18 

3.1 

12.9 

0.006 

(h) 

8.0 

2.93 

8.1 

15.7 

0.04 

Section  3.3.2.  Our  efficiency  decreases  rapidly  with  the  KBO  diameter:  from  nearly 
33%  at  77  =  30  km  to  2  x  10“®  at  77  =  0.5  km.  Note  that  at  77  =  30  km  we  ignore 
diffraction  effects  and  the  occultations  are  modelled  to  suppress  the  flux  completely 
for  several  consecutive  points,  depending  on  the  relative  velocity,  but  our  efficiency 
is  still  signihcantly  less  than  100%.  Some  of  our  lightcurves  are  too  noisy  to  allow 
detections. 

A  set  of  synthetic  events  recovered  by  our  pipeline  in  shown  in  Figure  3.7, 
and  the  parameters  of  each  plotted  event  are  given  in  Table  3.3.  Various  parameters 
affect  our  recovery  efficiency.  Our  efficiency  for  the  recovery  of  77  =  3.0  km  occulting 
KBOs  is  plotted  in  Figure  3.8  as  a  function  of  Mtaos;  SNR,  relative  velocity  Urei 
and  crowdedness  of  the  held. 

The  efficiency  as  a  function  of  SNR  is  plotted  in  Figure  3.8a  -a  few  targets 
at  SNR  >  100  are  left  out  of  the  plot.  The  efficiency  decreases  with  magnitude 
(Figure  3.8c)  by  about  a  factor  of  hve  between  magnitude  9  and  13.  The  dominant 
effect  here  is  the  decrease  in  SNR,  though  a  competing  effect  occurs  since  lower 
magnitudes  are  associated  with  larger  angular  sizes,  and  our  efficiency  decreases 
with  increasing  angular  size  (Figure  3.6).  Furthermore  there  are  many  more  dim 
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than  bright  stars  in  the  sky:  Figure  3.8f  shows  the  efficiency  as  a  function  of  Mtaos 
multiplied  by  the  number  of  TAOS  targets  at  that  magnitude.  The  highest  number 
of  detections  happen  for  stars  with  Mtaos  ~  12.5. 

The  behavior  of  our  efficiency  as  a  function  of  crowdedness,  where  the 
crowdedness  is  dehned  as  the  number  of  targets  in  that  TAOS  field  brighter  than 
-Rusno  =  13.5,  is  plotted  in  Figure  3.8b.  The  largest  number  of  detections  per  held 
is  achieved  for  helds  with  more  than  1500  Rusno  <  13.5  target  stars. 

Our  efficiency  as  a  function  of  the  relative  velocity  of  the  KBO  is  plotted 
in  Figure  3.8d.  Observing  at  a  pointing  where  the  relative  velocity  of  the  KBOs  is 
higher  boosts  the  event  rate  of  the  survey.  Our  efficiency  is  however  larger  for  smaller 
transiting  velocities,  particularly  for  small  KBOs  for  which  the  time-line  of  the  event 
is  shorter  than  one  of  our  data-points  at  opposition.  Ultimately,  the  effective  sky 
coverage  of  our  survey  depends  linearly  on  both  the  efficiency  and  the  velocity  (see 
section  4.2).  The  efficiency  multiplied  by  the  relative  velocity  Urei  is  plotted  against 
Urei  in  Figure  3.8f.  Pointing  near  opposition  increases  the  effective  coverage  of  our 
survey,  and  thus  it  increases  the  event  rate,  for  3  km  KBOs.  The  survey  strategy 
can  be  optimized  at  different  sizes  taking  into  account  the  size  dependent  efficiency 
X  Urei  as  well  as  the  expected  KBO  size  distribution. 
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Figure  3.8:  Recovery  efficiency  for  3  km  KBOs.  Panel  (a):  efficiency  as  function  of 
SNR.  Efficiency  versus  crowdedness  of  the  field  (b),  defined  as  the  number  of  targets 
brighter  than  RusNO  =  13.5.  Panel  (c):  efficiency  versus  magnitude  Mtaos  •  In  panel 
(e)  the  efficiency  as  a  function  of  magnitude  is  weighted  by  the  number  of  targets  at  that 
magnitude.  Panel  (d):  efficiency  versus  relative  velocity  of  the  KBO  targets.  In  panel  (f) 
the  efficiency  versus  relative  velocity  is  weighted  by  the  relative  velocity.  All  error  bars 
are  calculated  in  a  Poissonian  fashion  from  the  square  root  of  the  number  of  recoveries. 
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Chapter  4 


Constraints  on  models  of  the  Solar 
System  formation  and  evolution 
from  the  TAOS  data 


We  compute  the  number  of  events  expected  for  the  Kuiper  belt  formation 
and  evolution  models  by  Pan  &  Sari  (2005),  Kenyon  &  Bromley  (2004),  and  Be¬ 
navidez  &  Campo  Bagatin  (2009).  We  present  a  comparison  with  upper  limits  we 
derive  from  3.75  years  of  TAOS  data  which  constrains  the  parameter  space  of  these 
models.  This  is  the  hrst  direct  comparison  of  models  of  the  KBO  size  distribution 
with  data  from  an  occultation  survey.  Our  results  suggest  that  the  KBO  population 
is  indeed  strengthless  and  planetary  migration  played  a  role  in  the  shaping  of  the 
size  distribution^. 

^This  work,  together  with  Chapter  3,  was  accepted  for  publication  by  the  Astronomical  Journal 
under  the  title  The  TAOS  Project:  Upper  Bounds  on  the  Population  of  Small  KBOs  and  Tests 
of  Models  of  Formation  and  Evolution  of  the  Outer  Solar  System,  and  it  is  the  product  of  a 
collaboration  with  the  TAOS  team  and  with  M.  J.  Holman  (Harvard-Smithsonian  Center  for 
Astrophysics) 
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4.1  Introduction 


In  this  chapter  I  describe  the  constraints  on  the  KBO  population  derived 
from  our  analysis  of  3.75  years  of  TAOS  data,  discussed  in  Chapter  3.  In  Section  4.2 
I  derive  the  effective  coverage  of  our  survey  and  model-independent  limits  to  the 
number  of  objects  in  the  Kuiper  belt,  and  compare  these  results  with  those  of  similar 
surveys.  In  Section  4.3  I  compare  our  upper  limits  to  the  estimates  on  the  number 
of  KBOs  set  by  dynamical  simulations  for  JFC  progenitor  populations.  In  Section 
4.4  I  briefly  describe  models  for  the  formation  and  evolution  of  the  Kuiper  belts  and 
I  then  derive  and  discuss  constraints  to  these  models. 


4.2  Effective  coverage  and  upper  limits 


We  calculate  the  effective  sky  coverage  of  our  survey,  hie,  as: 


ne{D) 


1  Urel 

w{D)  ^  A  X 


(4.1) 


where  is  the  exposure  of  the  star  target  (the  duration  of  the  lightcurve  set), 
w{D)  the  weight  factor  for  that  diameter,  i.e.  the  fraction  of  lightcurves  implanted 
with  occultations  by  KBOs  of  diameter  D,  and  the  sum  is  carried  out  only  over 
the  lightcurves  where  events  are  recovered.  The  effective  coverage  of  our  survey, 
which  takes  into  account  our  efficiency,  is  plotted  in  Figure  4.1,  for  both  the  dataset 
published  in  Z08  (empty  squares)  and  for  the  current  work  (hlled  squares). 

We  can  use  Oe  to  calculate  model-independent  upper  limits:  at  each  size 
for  which  our  efficiency  calculation  was  conducted  we  calculate  the  number  density 
of  KBOs  from  the  TAOS  dataset  as  a  single-point  upper  limit.  These  limits  are 
shown  as  hlled  dots  in  Figure  4.2.  Each  point  represents  the  maximum  surface 
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Figure  4.1:  Effective  solid  angle  for  the  first  2  years  of  TAOS  data  (Z08,  empty  squares) 
and  for  the  current  3.75  year  dataset  (solid  line  and  filled  squares). 

density  of  KBOs  >  D.  Note  that  these  limits  are  not  statistically  independent  of 
each  other.  Figure  4.2  also  shows  the  similarly  calculated  upper  limits  of  Bickerton 
et  ah  (2008)  and  Bianco  et  al.  (2009a).  This  allows  a  direct  comparison  of  the  results 
from  our  survey  with  these  earlier  results  and  the  improvement  is  obvious:  nearly 
an  order  of  magnitude  at  700  m  and  over  an  order  of  magnitude  at  1  km  from  the 
most  recent  results,  those  of  Bianco  et  al.  (2009a). 

4.3  The  Jupiter  Family  Comets  progenitor  popu¬ 
lation 

The  Jupiter  Family  Comets  (JFCs)  are  believed  to  originate  in  the  Outer 
Solar  System.  In  this  scenario  the  giant  planets  generate  gravitational  perturbations 
that  affect  the  orbits  of  the  Outer  Solar  System  bodies,  injecting  them  into  the 
planetary  region,  where  they  are  captured  by  Jupiter.  The  orbital  inclination  of  the 
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Figure  4.2:  Upper  limits  from  the  TAOS  survey  (solid  dots).  Each  point  represents  the 
upper  limit  to  the  number  of  KBO  of  that  size  or  larger  given  the  TAOS  effective  coverage. 
Similarly  derived  upper  limits  from  Bickerton  et  al.  (2008)  at  1  km  (BKW,  empty  circles) 
and  Bianco  et  al.  (2009a)  at  1  km  and  0.7  km  (MMT,  empty  circles)  are  also  plotted. 
The  upper  limit  from  Bernstein  et  al.  (2004)  is  also  plotted  (HST).  The  estimates  of 
the  number  of  objects  in  the  Classical  Belt  (CB),  of  Plutinos  and  scattered  disk  objects 
(SD)  are  plotted,  as  derived  by  Levison  &  Duncan  (1997),  Morbidelli  (1997)  and  Volk  &: 
Malhotra  (2008)  respectively,  assuming  each  family  is  the  unique  precursor  of  JFCs. 
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JFCs  suggests  that  their  precursor  population  has  a  disk-like  distribution,  favoring 
thus  the  Kuiper  belt  over  the  Oort  Cloud  as  a  reservoir  (Volk  &  Malhotra,  2008, 
and  references  therein).  The  Classical  Belt  (CB),  the  Plutinos  and  the  scattered 
disk  (SD)  have  been  considered  as  precursors  in  various  studies.  The  dynamical 
characteristics  of  each  population  determines  the  efficiency  of  the  injection  process 
and  the  number  of  objects  in  each  progenitor  family  can  thus  be  derived  on  the  basis 
of  the  density  of  JFCs,  which  is  observationally  constrained  (see  Tancredi  et  ah  2006, 
and  references  therein).  Furthermore  the  size  distribution  of  JFCs  should  reflect  the 
size  distribution  of  the  progenitor  population.  Bernstein  et  ah  (2004)  and  Bianco 
et  ah  (2009a)  have  pointed  out  that  a  better  determination  of  the  size  distribution 
of  the  Kuiper  Belt  would  help  understand  the  origin  of  the  JFCs. 

We  assume  the  JFC  precursors  are  in  the  size  range  1  —  10  km  as  this 
is  observed  to  be  the  typical  size  of  JFCs  (Lowry  et  ah,  2008).  We  considered 
the  estimates  on  the  KBO  populations  (CB  and  Plutinos)  and  SD  derived  from 
dynamical  simulations  under  the  assumption  that  each  population  is  the  unique 
progenitor  of  the  JFCs,  and  we  compared  them  to  the  upper  limits  derived  from 
our  survey  (Figure  4.2).  We  use  the  estimate  of  Levison  &  Duncan  (1997)  for  a 
population  of  cometary  precursors  entirely  in  the  CB,  and  that  of  Morbidelli  (1997) 
for  Plutino  progenitors.  These  are  converted  into  a  surface  density  by  assuming 
for  each  population  a  projected  sky  area  of  10"^  deg^,  as  was  done  by  Bernstein 
et  ah  (2004).  We  consider  the  results  of  Volk  &  Malhotra  (2008)  for  a  progenitor 
population  in  the  SD  .  We  calculate  the  minimum  surface  density  of  SD  objects 
expected  in  the  region  of  sky  typically  observed  by  the  TAOS  survey.  For  this  we 
use  information  on  the  fraction  of  time  the  objects  spend  between  30  and  50  AU 
and  within  3°  of  the  ecliptic  plane  as  provided  by  Volk  &  Malhotra  (2008).  These 
estimates  on  the  number  of  objects  are  shown  in  Figure  4.2  as  horizontal  lines. 
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Our  results  rule  out  the  CB  as  unique  precursor  if  all  progenitors  have 
D  >  3  km.  Occultation  surveys  are  the  only  surveys  that  at  present  can  probe  this 
region  of  the  size  spectrum,  and  our  preliminary  result  shows  that  future  occultation 
surveys  will  be  able  to  derive  useful  constraints  on  the  origin  of  the  JFCs. 

4.4  Outer  Solar  System  collisional  models 

The  collisional  and  dynamical  evolution  of  the  Solar  System  shaped  the 
size  distribution  of  the  Kuiper  belt.  In  the  coagulation  scenario  the  Kuiper  belt 
was  hrst  shaped  out  of  collisions  among  the  primordial  planetesimals  (Kenyon  & 
Luu,  1999a).  The  belt  was  originally  populated  by  very  small  dust  grains,  with 
small  orbital  eccentricities  (e  <  0.01)  such  as  those  we  observe  in  circumstellar 
disks  around  other  stars  (Moro-Martin  et  ah,  2007).  Initially  these  small  objects 
merge  in  collisions  and  slowly  grow:  1  km  KBOs  in  the  Kuiper  belt  are  thus  formed. 
As  their  gravitational  cross  section  grows  larger  than  their  geometric  cross-section, 
gravitational  focusing  speeds  up  the  growth  rate  of  the  largest  bodies.  This  phase  is 
referred  to  as  runaway  growth,  and  objects  as  large  as  hundreds  of  kilometers  can 
form.  One  snch  population,  shaped  primarily  through  agglomeration  processes  is 
predicted  by  theory  to  have  a  power  law  distribution  in  diameter  dN/dD  oc 
with  power  gL  ~  4.5  (Kenyon  &  Lnn,  1999b).  Direct  observations  of  large  KBOs 
conhrm  the  power  law  behavior  in  this  regime,  the  gravitationally-dominated  region 
of  the  size  spectrum,  with  a  best  £t  of  gL  =  4.25  (Fraser  et  al.  2008,  Fuentes 
&  Holman  2008).  The  size  distribntion  of  these  large  objects,  for  which  gravity 
dominates  the  internal  strength,  is  remarkably  insensitive  to  parameters  snch  as 
Neptnne  stirring  or  the  internal  tensile  strength  of  the  KBOs. 

Meanwhile,  very  large  objects  in  the  planetary  region  of  the  Solar  Sys- 
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tern  are  also  forming  into  planets,  that  are  believed  to  undergo  significant  migra¬ 
tions  (Tsiganis  et  al.,  2005,  and  references  therein).  The  orbits  of  the  planetesimals 
are  then  stirred  up  via  gravitational  interaction  to  velocities  such  that  further  im¬ 
pacts  will  result  in  the  disruption  of  the  smaller  objects:  this  is  the  catastrophic 
collisions  phase  (Davis  &  Farinella,  1997;  Kenyon  &  Luu,  1999a;  Morbidelli  et  al., 
2008);  the  time  scale  to  reach  this  phase  is  estimated  to  be  between  10  Myr  and 
1  Gyr  (Kenyon  &  Bromley,  2001).  A  second  steady  state  occurs  for  very  small  ob¬ 
jects  (probably  tens  of  meters  and  smaller)  where  the  collisionally  evolved  population 
transitions  to  a  strength- dominated  regime.  Here  the  collisional  cascade  will  gener¬ 
ate  a  size  distribution  which  follows  a  power  law  with  index  gs  =  3.5  (Dohnanyi, 
1969;  Kenyon  &  Luu,  1999a),  also  in  a  fashion  that  is  largely  independent  on  the  de¬ 
tails  of  the  evolution  of  the  protoplanets.  Note  that  the  study  of  collisions  between 
icy  bodies  is  still  in  its  infancy,  and  future  work  in  this  held  will  permit  assessing 
the  behavior  of  colliding  small  strength-less  or  loosely  bound  particles  (Leinhardt, 
2008).  Similarly,  future  work  on  coupling  collisional  and  dynamical  evolution  codes, 
recently  pioneered  by  Charnoz  &  Morbidelli  (2007),  should  provide  further  insight 
in  the  behavior  of  the  size  distribution. 

We  will  refer  to  the  region  in  between  these  two  regimes  as  the  intermediate 
region.  The  extent  of,  and  the  size  distribution  in,  the  intermediate  region  are 
instead  very  sensitive  to  the  formation  and  evolution  parameters,  and  observational 
information  on  this  region  can  be  compared  to  evolution  models.  In  this  section, 
we  present  four  models  of  the  KBO  size  distribution  dN/dD,  three  models  from 
literature  and  a  simple  parametric  model.  On  the  basis  of  these  models  we  calculate 
the  number  of  events  expected  to  be  detected  by  the  TAOS  survey,  which  is  given 
by: 

A^exp  =  /  ^  dD. 
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Since  no  events  were  found  in  our  survey,  any  model  which  predicts  iVgxp  >  3  is  t 
ruled  out  by  TAOS  at  the  95%  conhdence  level. 

4.4.1  Pan  &  Sari  (2005) 

Pan  &  Sari  (2005,  hereinafter  PS05)  derived  a  fully  analytical  model  for 
the  size  distribution  of  KBOs  by  assuming  the  population  is  in  a  steady  state  and 
the  mass  is  constant  through  the  collisional  processes.  They  assume  for  most  of 
their  model  calculation  that  the  internal  strength  of  the  objects  is  negligible  (gravity 
dominated  objects).  This  assumption  is  motivated  by  studies  of  comets  and  asteroids 
(PS05  and  references  therein).  PS05  discuss  the  transition  to  the  fully  strength- 
dominated  regime,  where  they  model  the  size  distribution  as  a  power  law  with  power 
index  qs  =  3.5  according  to  Asphaug  &  Benz  (1996).  This  region  however  starts  at 
D  <  300  m  and  is  entirely  below  the  sensitivity  of  TAOS.  Under  these  assumptions 
PS05  derive  a  double  power  law  size  distribution  for  objects  D  >  300  m; 

dN/dD  ocd-‘^^  for  D  > 

dN/dD  ocd-®  for  D  <  D^.  (4.2) 

This  model  is  shown  in  Figure  4.3,  left.  The  slope  q  has  value  gr  =  5  for  large  objects 
and  qi  =  3  for  objects  in  the  intermediate  region.  They  are  thus  able  to  calculate 
self-consistently,  as  a  function  of  time,  the  location  of  the  break  the  power  D\^,  which 
represents  the  size  of  the  largest  KBOs  that  experienced  catastrophic  collisions.  The 
location  of  the  break  moves  toward  larger  objects  as  the  size  distribution  evolves; 
the  authors  predict  the  break  near  =  40  km.  Note  that  in  the  intermediate 
region  between  D]^  and  the  second  break,  which  is  the  region  that  TAOS  can  probe,  a 
realistic  size  distribution  is  expected  to  have  an  oscillatory  behavior  while  preserving 
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Figure  4.3:  PS05  model  and  expected  event  rate  for  TAOS.  Left.  Triangles  are  the  data 
from  Fraser  &:  Kavelaars  (2008).  Empty  circles  are  the  data  from  Bernstein  et  al.  (2004). 
We  use  these  direct  observation  data  to  set  the  location  of  the  large  end  size  distribution. 
The  model  is  parametrized  with  two  slopes,  ql  =  5  and  qi  =  3.  The  positions  of  the  break 
are  at  T>b  =100,  80,  40,  20,  and  10  km.  Right.  Expected  yield  of  events  in  3.75  years  of 
TAOS  data.  The  horizontal  dashed  line  represents  the  highest  number  of  events  allowed 
by  no-detections  (3  events  at  the  95%  confidence  level).  Any  model  above  this  limit  is 
ruled  out  by  our  survey. 

an  average  slope  of  qj  =  —3  (PS05). 

We  parametrized  the  PS05  model  with  two  slopes  and  a  hard  junction 
point.  Distributions  with  break  points  Db  =  100,  80,  40,  20,  and  10  km  are 
shown  in  Figure  4.3.  The  left  panel  shows  our  parametrization  of  the  model,  and 
the  corresponding  predicted  number  of  events  for  the  TAOS  data  analyzed  in  this 
paper  is  plotted  in  the  right  paneh. 

On  the  basis  of  our  no-detections  result,  break  diameters  of  Db  <51.3  km 
are  excluded  at  95%  conhdence  level.  This  is  consistent  with  the  data  from  direct 
observations  and  with  the  authors’  interpretation  of  the  model:  in  absence  of  stirring 

^Note  that  the  uppermost  curve  in  the  right  panel  of  Figure  4.3,  which  corresponds  to  a  simple 
extension  of  the  large-size  power  law  distribution,  represents  the  naive  expectation  before  Bernstein 
et  al.  (2004)  appeared.  This  was  the  initial  design  target  of  the  TAOS  project.  TAOS  has  met  and 
exceeded  its  original  design  goals! 
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Figure  4.4:  Our  modeling  of  the  KB04  results.  The  range  of  results  of  the  simulations  of 
KB04  is  represented  by  the  shaded  region.  A  few  of  our  models  are  plotted  as  solid  lines. 
All  models  are  normalized  by  the  slope  of  the  large  KBO  size  distribution  4.5. 

by  Neptune,  the  location  of  the  break  suggests  that  the  KBO  population  is  comprised 
of  objects  with  little  internal  strength. 


4.4.2  Kenyon  &:  Bromley  (2004) 

Bernstein  et  al.  (2004):  here  collisions  destroy  weak  KBOs  and  models 
with  Neptune  stirring,  or  weakly  bound  KBOs  produce  a  more  significant  dip.  This 
feature  is  followed  by  an  excess  with  respect  to  the  nominal  power  law  for  2-15  km 
KBOs  which  varies  in  amplitude  between  a  factor  of  a  few  and  a  factor  of  a  few  tens, 
depending  on  the  internal  strength  and  presence  of  Neptune  (see  Figures  4.4  and 
Figures  4.5  &  4.6,  left  column).  The  size  distribution  remains  sensitive  to  the  details 
of  the  models  down  to  about  H  =  50  m,  where  once  again  a  power  law  behavior 
begins,  with  power  3  <  g  <  5  in  the  strength-dominated  regime. 

We  parametrize  the  KB04  suite  of  models  to  allow  us  to  compare  them 
directly  to  our  survey.  Specihcally,  we  model  the  large  size  distribution  as  a  power 
law  with  gL  =  4.5  and  the  intermediate  region  with  gi  =  2.  We  also  model  the 


excess  {bump)  around  D  =  5  km  with  a  Gaussian,  so  that: 
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(4.3) 


We  hx  the  width  of  the  Gaussian  excess  to  Ue  =  3.5  km,  we  set  the  location  to 
/iex  =  5.5  km  or  pex  =  1-6  km.  The  intensity  of  the  excess  is  determined  by  /ex;  we 
consider  models  with  an  excess  of  /ex  =  0,  10  and  100.  Note  that  as  the  break 
diameter  moves  towards  large  sizes  the  models  naturally  simulate  the  small  dip  near 
20  km.  Our  results  are  not  sensitive  to  the  presence  of  this  depletion. 

In  Figure  4.4,  the  range  of  results  of  the  KB04  simulations  is  shown  as  a 
shaded  region  and  a  few  parameter  choices  for  our  model  are  shown  as  solid  lines. 
The  size  distribution  is  shown  scaled  by  the  slope  of  the  large  size  region  (/1“^'^).  Fo¬ 
cusing  on  the  region  near  the  transition  between  primordial  and  collisionally  evolved 
population,  an  excess  near  5  km  and  a  depletion  near  25  km  are  both  visible  and  well 
represented  in  our  models.  The  behavior  of  the  simulations  at  D  >  100  km  is  not 
well  £t  by  our  models  but  this  does  not  affect  our  event  rate  calculation,  since  the 
expected  number  of  occultation  by  such  large  objects  is  negligibly  small.  Figure  4.5 
shows  the  models  with  pex  =  5.5  km  (left)  and  the  corresponding  expected  number 
of  events  for  the  TAOS  dataset  (right)  in  absence  of  a  bump  (top)  and  for  bump 
intensities  Ax  =  10  (center).  Ax  =  100  (bottom).  Figure  4.6  shows  the  models  with 
Pex  =  1-6  km  (left)  and  the  corresponding  event  rates  for  the  TAOS  dataset  (right) 
for  a  bump  intensities  Ax  =  10  (top)  and  100  (bottom).  Both  figures  show  models 
with  break  locations  D\^  =  60,  40,  20,  10,  and  with  no  break. 

For  a  bump  centered  at  pex  =  5.5  km  break  points  /lb  >  75.3  km  are 
allowed  for  any  combination  of  parameters.  We  can  constrain  the  intensity  of  the 
excess:  any  break  diameter  smaller  that  75.3  km  is  ruled  out  for  models  with  a 
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/ex  =  100,  break  diameters  -Db  <  32.6  km  are  ruled  out  for  /ex  =  10  and  breaks 
at  /lb  <  16.5  km  are  ruled  out  in  absence  of  the  excess.  When  moving  the  break 
towards  smaller  sizes  (/Xex  =  1-6  km),  break  locations  are  ruled  out  for  /lb  <  63.3  km 
if  /ex  =  100,  and  /lb  <  28.3  km  if  /gx  =  10.  Note  that  the  KB04  simulations  show 
that  models  with  weaker  KBOs  and  Neptune  stirring  favor  a  location  of  the  bump 
at  smaller  sizes. 

All  together  our  result  strongly  favors  models  that  incorporate  the  effects 
of  Neptune  stirring  and  weaker  KBOs,  where  the  bulk  strength  Qb  is  roughly  Qb  < 

103  grg  g-i 

4.4.3  Benavidez  &:  Campo  Bagatin  (2009) 

Recent  simulations  by  Benavidez  &  Campo  Bagatin  (2009,  BCB09  here¬ 
inafter  )  divide  the  Kuiper  belt  into  three  dynamical  families  -  the  CB,  the  Plutinos 
(Resonant  Population)  and  the  SD  -  and  follow  the  collisional  evolution  of  each, 
while  taking  into  account  the  physics  of  the  fragmentation  of  icy  and  rocky  bodies 
at  the  typical  relative  velocities  of  KBOs.  This  suite  of  models  ignores  the  effects 
of  Neptune  stirring.  The  models  incorporate  four  scaling  laws  for  fragmentation:  a 
simple  scaling  driven  by  gravitational  self-compression,  two  scaling  laws  which  in¬ 
clude  both  self-compression  and  the  effects  of  strain-rate,  as  described  by  Farinella 
et  ah  (1982),  with  different  diameter  dependency  and  and  one  that 

follows  the  modeling  of  Benz  &  Asphaug  (1999)  for  icy  bodies.  They  then  vary  the 
material  strength  of  the  objects  between  10^  and  lO^erg  cm“^. 

They  predict  a  power  law  size  distribution  for  large  objects,  with  slope  gL  ~ 
4.0,  and  with  a  first  break  around  100  km.  The  size  distribution  then  departs  from 
power  law  behavior  for  all  parameter  choices,  with  two  to  four  orders  of  magnitudes 
fewer  KBOs  in  the  ~  1  km  region  than  the  nominal  power  law  would  predict.  The 


size  distribution  then  follows  again  a  power  law,  with  slope  qs  =  3.5,  in  the  strength- 
dominated  regime.  The  BCB09  models  are  shown  in  Figure  4.7,  left.  We  simplihed 
the  size  distribution  behavior  in  the  intermediate  region  with  a  single  power  law, 
though  typically  the  behavior  is  oscillatory. 

These  models  are  all  allowed  by  the  TAOS  data,  partly  because  of  the  loca¬ 
tion  of  the  initial  break  at  =  100  and  the  slope  of  the  large  end  size  distribution. 
The  break  location  is  at  the  large  end  of  what  is  allowed  by  direct  observations 
(indeed  outside  of  the  Fraser  &  Kavelaars  (2008)  allowed  range  of  Hb  G  [50,  95  km] 
with  assumption  of  a  6%  albedo).  The  location  of  this  first  break  does  not  evolve  in 
the  BCB09  simulations  from  what  is  set  as  an  initial  condition.  Similarly  the  slope 
at  the  large  end  of  the  size  distribution  gL  =  4.0,  which  is  slightly  shallower  than  the 
current  best  £t  value  from  Fraser  &  Kavelaars  (2008)  and  Fuentes  &  Holman  (2008) 
is  input  as  an  initial  condition  and  the  value  hardly  changes  in  the  simulation.  As  the 
parameters  relative  to  the  large  size  end  of  the  distribution  are  more  firmly  pinned 
down  by  direct  observations,  future  occultation  data  will  place  stronger  constraints 
on  the  details  of  the  shape  of  the  size  distribution  below  the  first  break  and  thus 
the  details  of  the  fragmentation  mechanisms,  providing  information  on  the  internal 
structure  of  the  KBOs. 

4.4.4  Generic  3— regime  model:  constraints  on  the  interme¬ 
diate  region  of  the  size  spectrum 

While  direct  surveys  can  set  constraints  on  the  slope  of  the  size  distribution 
for  large  KBOs  (gb)  and  location  of  the  first  break  in  the  power  law  (T*bi),  our 
survey  provides  the  only  current  dataset  able  to  specifically  probe  what  happens 
after  this  first  break.  A  simple  3-regime  model  allows  us  to  parametrize  separately 


the  primordial  region,  the  intermediate  region  and  the  fully  collisional,  strength- 
dominated  regime.  Knowing  that  both  the  large  (gravitationally-dominated)  and 
the  small  (strength-dominated)  regions  of  the  size  spectrum  are  rather  insensitive 
to  the  details  of  the  internal  structure  and  evolution  of  the  Kuiper  belt,  it  is  the 
intermediate  region  that  contains  the  most  information  about  the  physical  details 
of  the  KBOs.  We  grossly  simplify  the  expected  structure  in  the  intermediate  regime 
and  describe  it  with  a  power  law.  The  parameters  of  one  such  models  would  then 
be  three  slopes  gL,  gi,  and  gs,  and  two  break  locations  Dbi  and  D\,2. 

We  set  the  slope  of  the  large  end  of  the  size  spectrum  and  the  location  of 
the  first  break  to  the  best  £t  to  the  data  from  direct  observations:  gL  =  4.25  and 
Dhi  =  75  km.  We  model  the  intermediate  region  as  a  plateau,  or  with  a  shallow 
slope:  gi  G  [0,  3],  and  the  small  size  end  of  the  spectrum  as  a  power  law  with  slope 
3.5  for  the  strength-dominated,  collisionally  evolved  population.  Figure  4.8  shows 
our  3-regime  model  (left)  and  the  limit  we  can  set  to  the  intermediate  slope,  gi,  and 
second  break  location,  Zi)b2,  phase  space  (right).  Any  pair  of  values  gi  and  D\^2  that 
fall  above  the  solid  line  are  excluded.  At  this  time  we  can  not  only  use  this  result 
as  a  guideline  in  survey  design  but  also  for  the  comparison  of  any  size  distribution 
model  with  the  present  TAOS  data,  particularly  as  the  data  from  direct  surveys 
help  constrain  more  firmly  the  large  end  size  slope  and  location  of  the  first  break. 


90 


R  magnitude 
25  30 


R  magnitude 


R  magnitude 
25  30 


R  magnitude 
25  30 


R  magnitude 


R  magnitude 


Figure  4.5:  KB04  models.  The  left  column  shows  the  differential  size  distribution, 
parametrized  as  two  slopes  (^l  =  4.5  and  qi  =  2)  and  on  the  right  side  are  the  corre¬ 
sponding  event  rates  for  the  TAOS  survey.  Symbols  are  the  same  as  in  Figure  4.3.  The 
top  row  shows  the  models  without  bump,  the  center  and  bottom  panel  with  a  Gaussian 
bump  centered  at  ^ex  =  5.5  km  of  intensity  /ex  =  10  (top),  and  /ex  =  100  (bottom),  with 
/lb  =60,  40,  20  and  10  km,  and  no  break. 
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Figure  4.6:  Kenyon  and  Bromley  model,  as  in  Figure  4.5  but  with  Gaussian  bump  centered 
at  Hex  =  1-6  km,  Igx  =  10,  top,  and  lex  =  100,  bottom.  Right:  corresponding  number  of 
events  in  our  TAOS  survey.  Symbols  are  the  same  as  in  Figure  4.3 
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Figure  4.7:  Our  parametrization  of  BCB09  models:  the  models  differ  in  the  prescription 
for  fragmentation  and  all  models  are  parametrized  as  a  series  of  three  slopes.  Symbols 
are  the  same  as  in  Figure  4.3.  The  first  break  points  is  fixed  at  100  km.  The  slope  for 
the  smallest  size  objects  is  set  to  3.7.  The  location  of  he  second  break  and  intermediate 
slope  are  5  km  and  1.0,  3.6  km  and  2.0,  4.6  km  and  2.5  and  0.36  km  and  3.0.  Right: 
corresponding  number  of  events  in  our  TAOS  survey. 
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Figure  4.8:  Three-slope  model  described  in  Section  4.4.4  (left)  for  the  minimum  second 
break  location  i4b2  allowed  by  slopes  qi  =  0  and  3.  Slope-break  location  phase  space  for 
the  three-slopes  model  (right).  The  region  above  the  line  is  excluded  by  our  data. 
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Chapter  5 


Exploring  the  Solar  System 
beyond  the  Kniper  belt 


We  explored  5  x  10^  star  hours  of  data  from  the  Taiwanese  American 
Occultation  Survey  (TAOS),  the  dataset  described  in  Chapter  3,  searching  for  oc- 
cultations  of  bright  stars  by  solar  system  bodies  in  orbit  between  100  and  1000  AU. 
We  found  no  such  event.  Our  hnal  goal  will  be  to  compare  our  no-detection  result 
with  predictions  from  models  of  the  formation  of  the  Sedna  population  and  of  the 
scattered  disk,  thus  setting  upper  limits  on  the  number  density  of  Sedna-like  and 
scattered  disk  objects  and  constraints  on  formation  evolution  scenarios^. 

5.1  Introduction 

The  outer  solar  system  is  populated  by  icy  bodies  that  formed  in  the  early 
stages  of  planetary  formation.  While  direct  imaging  surveys  can  probe  the  classical 
belt  and  the  resonant  population  of  the  Kuiper  belt  (KB)  down  to  objects  as  small 
as  tens  of  kilometers,  the  outer  solar  system  extends  well  beyond  50  AU,  in  various 
^The  work  presented  in  this  chapter  is  the  product  of  a  collaboration  with  the  TAOS  team. 
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dynamical  families  of  objects:  the  scattered  disk  (SD),  extended  scatter  disk  (ESD) 
and  inner  and  outer  Oort  cloud.  In  these  regions  objects  as  large  as  tens,  hundreds 
or  even  thousands  of  kilometers  in  diameter  would  elude  direct  detection. 

With  the  term  SD  we  generally  refer  to  the  region  of  orbital  space  with 
perihelion  g  >  30  AU  and  semi-major  axis  a  >  50  AU  (see  for  example  Gomes  et  ah 
2005b).  The  inclination  distribution  of  SD  objects  is  larger  than  that  of  the  classical 
KB,  but  generally  limited  to  i  <  30°.  These  objects  originate  in  the  classical  KB: 
with  pericenter  just  outside  the  orbit  of  Neptune  their  orbits  evolved  on  timescales  of 
Gyr  through  repeated  long  distance  interactions  with  Neptune  (Duncan  &  Levison, 
1997). 

With  larger  inclination  and  eccentricity  distribution  than  the  SD,  and  ex¬ 
tending  farther  out  in  the  outer  solar  system,  the  origin  of  the  ESD  is  not  certain. 
We  generally  classify  object  as  ESD  objects  if  g  >  38  AU,  thus  decoupling  them 
from  Neptune.  Mechanisms  involving  mean- motion  trapping  in  conjunction  with 
Kozai  chaos  have  been  proposed  (Gomes  et  ah,  2005b;  Volk  &  Malhotra,  2009),  as 
well  as  early  stellar  encounters,  to  explain  the  evolution  of  planetesimals  from  the 
KB  or  SD  into  the  ESD. 

The  long  period  comets  that  cross  the  inner  solar  system  provide  observa¬ 
tional  evidence  of  the  existence  of  an  isotropic  family  of  objects  extending  farther 
out  in  the  solar  system.  The  Oort  cloud  formed  as  gravitational  interactions  with 
the  giant  planets  pushed  the  aphelia  of  planetesimals  from  the  SD  into  the  region  of 
influence  of  the  Galactic  tide  (a  >  3,  000  AU,  Duncan  et  al.  1987).  The  Oort  cloud 
comprises  two  distinct  structures:  the  inner  Oort  cloud,  extending  to  20,  000  AU, 
and  the  outer  Oort  cloud,  a  spherical  structure  extending  as  far  out  as  tens  of  thou¬ 
sands  of  AU.  The  long  period  comets  originate  from  the  latter,  when  the  orbit  of 
and  Oort  cloud  object,  being  far  from  the  influence  of  the  Sun,  is  perturbed  by  star 
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passages.  The  outer  Oort  cloud  is  loosely  bound  to  the  Sun,  thus  a  signihcant  loss 
rate  is  to  be  expected;  the  inner  Oort  cloud  should  replenish  the  comets  lost  from 
the  outer  cloud,  though  mechanisms  for  this  process  are  uncertain  and  might  require 
interactions  with  external  bodies  (Brasser  et  ah,  2006). 

The  discovery  of  Sedna  (2003  VB12,  Brown  et  al.  2004)  suggested  the  in¬ 
triguing  possibility  of  the  existence  of  a  distant  population  of  objects  with  perihelia 
as  close  as  a  few  hundred  AU  from  the  Sun,  but  that  spend  most  of  their  time  far  from 
the  region  of  dynamical  influence  of  Neptune.  Sedna  ha  s  perihelion  q  =  76.156  AU, 
semi-major  axis  a  =  525.86  AU,  eccentricity  e  =  0.855  and  inclination  i  =  11.934°. 
No  mechanism  yet  has  been  found  to  place  Sedna  on  its  orbit  without  invoking  the 
presence  of  an  unobserved  perturbing  object:  an  unknown  rogue  planet  or  a  pass¬ 
ing  star.  Continuing  direct  and  indirect  surveys  have  detected  no  other  object  like 
Sedna  (Schwamb  et  al.  2009,  and  references  therein,  Wang  et  al.  2009  -  hereinafter 
W09). 

Morbidelli  &  Levison  (2004)  have  excluded  scenarios  where  Sedna  evolves 
to  its  current  orbital  parameters  from  the  SD  or  KB  under  the  tidal  influence  of 
a  massive  early  cold  trans-Neptunian  disk.  Similarly  unlikely  are  scenarios  where 
Sedna  evolves  under  the  influence  of  Neptune  in  an  early  phase  when  Neptune’s 
orbit  had  a  higher  eccentricity.  This  scenario  acts  on  time  scales  that  are  long, 
while  Neptune  today  has  a  nearly  circular  orbit.  Also  this  mechanism  generates  a 
larger  nnmber  of  objects  in  the  50  <  a  <  100  AU,  than  a  >  100  AU,  while  both 
the  detached  object  that  we  know  populates  the  region  a  >  200  AU  (Sedna  and 
2000  CR105).  Gladman  &  Chan  (2006)  hnd  that  the  presence  of  rogue  planets  in  the 
trans-Neptnnian  region,  ejected  from  the  solar  system  on  time  scales  of  ~  200  Myr, 
could  raise  the  perihelia  of  the  SD,  thus  generating  an  ESD,  and  it  could  also  produce 
objects  with  large  orbital  inclinations  like  2004  XR190,  2003  UB313,  2000  CR105, 


96 


and  Sedna. 


Another  promising  scenario  sees  the  solar  system  forming  in  an  embed¬ 
ded  star  cluster.  The  presence  of  passing  stars  may  explain  at  once  the  orbit  of 
Sedna  (Morbidelli  &  Levison,  2004;  Brasser  et  ah,  2007;  Kaib  &  Quinn,  2008),  the 
replenishing  rate  of  the  outer  Oort  cloud  and  the  truncation  of  the  Classical  KB  at 
50  AU  (Brasser  et  ah,  2007,  and  references  therein). 

Each  one  of  these  formation  and  evolution  scenarios  leaves  a  distinct  mark 
in  the  dynamical  characteristics  of  the  outer  solar  system  families.  Disentangling 
the  different  scenario  requires  however  that  we  further  our  observational  knowledge. 
Probing  these  remote  regions  of  the  solar  system  may  profoundly  impact  our  un¬ 
derstanding  of  solar  system  formation,  revealing  not  only  details  of  the  interaction 
among  early  solar  system  objects,  but  as  well  details  of  the  Sun’s  birth  environment. 

Occultation  data  from  the  TAOS,  Taiwanese  American  Occultation  Survey, 
have  been  used  to  set  upper  limits  to  the  population  of  objects  in  Sedna-like  orbits 
(W09).  TAOS  was  designed  to  probe  the  Kuiper  Belt  population  of  objects  smaller 
than  the  current  observational  limit  (diameter  D  <  10  km).  With  simultaneous 
multi-telescope  high  cadence  observations  the  data  collected  by  TAOS  allowed  us  to 
produce  the  strongest  upper  limits  on  the  small  end  of  the  KBO  population  (Zhang 
et  ah,  2008;  Bianco  et  ah,  2009b). 

Here  we  present  an  analysis  in  progress  of  5  x  10®  star  hours  of  three- 
telescope  4  or  5  Hz  observations  collected  between  January  2005  and  August  2009, 
data  described  in  Chapter  3.  These  data,  a  substantially  larger  dataset  than  the 
one  analyzed  in  W09,  are  being  scanned  to  detect  shallow  flux  drops  affecting  a  few 
to  a  few  tens  consecutive  data  points  simultaneously  in  all  of  our  telescopes. 
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5.2  Search  algorithm 


The  signature  of  an  occultation  of  a  bright  star  by  a  body  at  A  =  100  — 
1000  AU,  would  appear  in  a  TAOS  lightcurve  as  a  shallow  flux  drop  affecting  5-50 
consecutive  points,  depending  on  the  relative  velocity  and  size  of  the  occulter,  and 
possibly  displaying  prominent  diffraction  features.  We  want  to  scan  our  lightcurves 
looking  for  statistically  significant  deviations  of  the  flux  from  the  local  baseline. 
We  define  two  windows:  a  larger  baseline-window  W  over  which  calculate  our  local 
mean,  /^y,  and  local  standard  deviation,  and  a  smaller  event-window  w,  of  size 
comparable  to  the  duration  of  the  events  we  are  searching  for.  We  select  W  =  61 
time-stamps  and  w  =  11,  17,  and  21  time-stamps. 

We  hrst  hlter  our  lightcurves,  subtracting  from  each  point,  i,  the  local 
mean  fw  of  the  surrounding  region,  weighted  by  local  the  standard  deviation  aw- 


fw  {i)  -  X]  /i  +  X] 


W 


i=*+f 


aw  = 


h{i)  = 


W-w 


^  {fj  -  fw)  +  ^  {fj  -  fw)  , 


f{i)  -  fw{i). 


awW 


(5.1) 

(5.2) 

(5.3) 


h{i)  represents  the  deviation  from  the  mean  in  a  baseline-window  W  surrounding  a 
potential  event,  but  disregarding  the  central  w  time-stamps  that  would  be  impacted 
by  the  event  itself.  The  deviation  from  the  mean  is  here  expressed  in  units  of 
standard  deviation  aw,  where  aw  is  also  calculated  over  the  baseline-window  W, 
disregarding  the  central  w  points. 

For  each  event-window  size  we  rank  the  new  lightcurves  h{i)  for  each  tele- 


scope  (A,  B  and  D,  see  Chapter  3)  from  the  lowest  to  the  highest  value,  so  that 
the  time-stamp  with  largest  deviation  from  the  mean,  where  the  flux  is  suppressed, 
would  be  ranked  r*  =  1.  This  algorithm  is  similar  to  that  used  in  W09,  but  it 
accounts  for  the  non-stationary  nature  of  the  noise  in  our  data.  The  algorithm  we 
implemented  here  thus  allows  us  to  calculate  exactly  the  p-value  of  each  of  our  rank 
triplet:  set  of  three  ranks  associated  with  the  same  time-stamp  for  our  three  tele¬ 
scopes.  In  absence  of  events,  just  like  we  did  in  Chapter  3.3.2,  we  can  calculate  the 
probability  distribution  of  the  quantity: 


where  A^t  is  the  number  of  telescopes  simultaneously  observing  the  star,  Np  the 
number  of  points  in  the  lightcurve,  and  rik  the  rank  of  the  f-th  point  in  the  lightcurve 
observed  by  telescope  k,  and  we  can  compute  the  probability  of  a  random  variable 
Z  arising  from  this  distribution:  P{Z  >  z{i))  =  Given  the  number  of  points  in 
our  dataset  we  expect  <  0.2  false  positive  events  generated  by  noise  fluctuations  for 
a  threshold  =  3.0  x  10“^^  (see  Section  3.3.2). 

Figure  5.1,  top,  shows  occultation  lightcurves  of  objects  of  D  =  10  km 
(left)  and  D  =  5  km  (right)  diameter  at  A  =  400  AU,  occulting  magnitude  i?  ~  12 
stars.  The  bottom  panel  shows  the  quantity  z{i)  for  true  star  lightcurves  implanted 
with  the  synthetic  events  shown  above,  for  an  event-window  size  w  =  10.  An  event 
generated  by  an  object  as  small  as  5  or  10  km  is  obviously  visible  at  A  =  400  AU. 
These  objects  are  far  beyond  the  limits  of  direct  observations! 


99 


1 


Figure  5.1:  Top  panel:  occultation  signature  of  a  H  =  10  km  object  (left)  and  a  D  =  5  km 
(right)  object  at  A  =  400  AU  occulting  a  magnitude  i?  ~  12  star.  The  impact  parameter 
is  assumed  to  be  6  =  0.  Bottom  panel:  product  of  the  ranks  of  the  lightcurves  of 

telescopes  A,  B,  and  D  in  the  region  implanted  with  the  occultation. 

5.3  Event  rate  calculation  for  Sedna-like  objects 


We  are  interested  in  assessing  the  density  of  objects  in  Sedna-like  orbits. 
We  will  test  our  recovery  efficiency  for  objects  in  the  size  range  1  <  H  <  30  km, 
at  distances  of  100,  200,  300,  500  and  1000  AU. 

For  every  diameter  D  and  distance  A  in  which  we  are  interested,  we  want 
to  calculate  the  expected  event  rate  for  our  survey.  This  event  rate  will  depend 
on  the  relative  velocity  nrei(e,A)  of  the  occulter,  where  e  is  the  elongation  of  the 
occulter.  To  hrst  approximation  use  Equation  1.8  to  estimate  the  relative  velocity 
of  the  occulter.  This  formula  neglects  the  proper  motion  of  the  object.  We  plan  on 
eventually  including  specihc  models  for  the  orbital  distribution  of  both  Sedna-like 
and  SD  objects,  thus  introducing  specihc  velocity  distributions  in  our  estimate  of 

Uel- 

For  a  given  star  lightcurve  (*),  and  an  observing  angle  e*,  we  calculate  the 
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expected  event  rate  i?*  of  occultations  as: 


d^R,{D,  A)  dMD,  A)  t;,ei(e,,  A)H,{D,  A) 
dD  dA  ~  dD  dA  A2 


where  H^{D,  A)  is  the  event  cross  section,  set  to  the  diameter  of  the  first  Airy  ring 
(see  Nihei  et  ah  2007  and  Section  1.3.1  of  this  work),  and  where  n{D,A)  is  the 
number  density  of  objects  with  diameter  D  at  distance  A,  which  is  what  we  want 
to  constrain  based  on  the  number  of  events  detected  by  our  survey. 

The  number  of  events  expected  for  our  survey  is 


d^N,^^{D,A) 

dD  dA 


E 


d^R,{D,A) 

dD  dA 


E^ri^{D,A), 


(5.6) 


with  the  exposure  (the  duration  of  the  data  run)  and  ri^{D,A)  is  the  efficiency 
of  recovery  of  objects  of  diameter  D  at  distance  A  in  the  lightcurve  we  are  analyz¬ 
ing,  and  the  sum  is  conducted  over  all  our  lightcurves.  Since  we  implant  exactly 
one  occultation  -  by  an  occulter  of  diameter  D  at  distance  A  -  in  each  one  of 
our  lightcurves,  the  efficiency  of  recovery  assumes  value  of  =  1  if  the  event  is 
recovered,  or  ?]*  =  0  if  it  is  not.  Denoting  the  fraction  of  lightcurves  implanted 
with  objects  of  diameter  D  with  wd,  and  with  objects  at  distance  A  with  wa, 
Equation  5.6  takes  the  form: 


d^Aexp(D,A) 

dD  dA 


d^n(D,A)  1 

dD  dA  wdWa 


E 

*,rec 


Vr:el{e^,A)H^{D,A)E^ 

A2 


(5.7) 


where  this  time  the  sum  is  conducted  only  over  the  lightcurves  with  recovered  events. 
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Figure  5.2:  Minimum  coverage  of  our  survey:  the  coverage  of  W09  is  projected  to  account 
for  the  increased  exposure  only. 

As  seen  in  Section  4.2  of  this  work  then: 


1  A)i?* 

wdWa  ^  A^ 

* 

is  the  effective  coverage  of  our  survey. 

Compared  to  the  W09  the  exposure  of  our  survey  has  increased  by  a  factor 
2.95.  From  preliminary  tests  we  expect  the  modihed  search  algorithm  we  use  here 
to  increase  the  efficiency  of  our  search,  but  we  are  unable  to  precisely  quantify  this 
improvement  at  this  time.  Figure  5.2  shows  our  effective  coverage  as  a  function 
of  diameter  D  separately  for  all  A  distances  tested  accounting  for  the  increased 
exposure  but  ignoring  the  increase  in  efficiency. 


(5.8) 


f^e(/^,A) 
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5.4  Renewed  limits  to  the  Sedna  population 


In  W09  we  set  upper  limits  to  the  population  of  Sedna-like  objects.  These 
limits  are  reproduced  here  accounting  for  the  increased  coverage  of  our  present 
survey,  but  disregarding  the  expected  increase  in  efficiency.  These  are  therefore 
conservative  upper  limits  to  the  Sedna  population. 

We  derive  constraints  to  the  number  density  of  Sedna-like  objects  sepa¬ 
rately  for  every  distance  A  we  tested.  We  have  to  assume  a  functional  form  for 
the  diameter  dependence:  a  size  distribution.  We  assume  a  simple  power  law  size 
distribution  n  oc  D^~'^  so  that: 


rj  79 

—  =  no{q  -  (5.9) 

with  uq  the  number  of  objects  of  diameter  D  >  1  km  at  distance  A.  We  expect  in 
our  survey 


rDmax 

Ae,,p(A)  =  no  (g  -  1)  /  A)  dD,  (5.10) 

^  -^min 

events  from  objects  at  distance  A.  Because  we  recovered  no  events  we  can  exclude 
at  95%  conhdence  level  any  value  of  q  that  predicts  more  then  3  events.  For  values 
of  q  ranging  from  q  =  3.0  to  g  =  6.0  we  thus  set  upper  limits  on  the  number  density 
of  objects  at  distance  A  in  Figure  5.3:  models  above  the  plotted  lines  are  excluded. 

In  WOO  we  derived  a  lower  limit  to  the  density  of  Sedna-like  objects,  based 
on  the  existence  of  Sedna:  any  value  of  g  that  predicts  fewer  than  0.05  objects 
larger  than  Sedna  is  excluded  at  the  95%  conhdence  level.  This  correspond  to  a 
minimum  number  density  of  Sednas  (i.e.  objects  larger  than  D  =  1600  km  at 
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Figure  5.3:  Maximum  number  density  of  D  >  1  km  objects  as  a  function  of  the  slope  of 
a  power  law  size  distribution.  Distances  of  100,  200,  300,  500,  and  1000  AU  are  assumed. 
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Figure  5.4:  Maximum  number  of  density  (left  axis)  and  absolute  number  (right  axis)  of 
objects  larger  than  Sedna  {D  >  1600  km)  as  a  function  of  the  slope  of  a  power  law  size 
distribution  (solid  line).  The  region  above  the  line  is  excluded  by  the  TAOS  result.  The 
lower  limit  set  by  the  very  existence  of  Sedna  is  also  plotted  (dashed  line) 

100  AU)  Us  <  1.2  X  10“®  deg~^.  Combining  this  lower  limit  with  our  upper  limits 
we  can  state  that  slopes  q  >  5.2  are  excluded  at  95%  c.l.  (see  Figure  5.4). 


5.5  Future  work 

So  far  we  reproduced  the  work  done  in  W09,  taking  advantage  of  our 
increased  coverage.  We  expect  to  be  able  to  derive  more  stringent  limits  due  to  the 
increased  efficiency  of  our  search  algorithm.  More  interestingly,  though,  we  want 
to  include  possible  velocity  and  orbital  distributions  for  both  the  Sedna-like  and 
Scattered  Disk  objects. 


105 


Chapter  6 


Conclusions 


I  collected  and  analyzed  a  large  amount  of  data  looking  for  serendipitous 
occupations  of  bright  stars  by  objects  in  the  outer  solar  system.  Occultations  would 
allow  us  to  detect  the  presence  of  objects  too  small  or  too  distant  to  be  observed 
directly.  I  found  no  occultations  by  outer  solar  system  objects  in  my  datasets  and  I 
am  able  to  set  upper  limits  on  the  KBO  population  and  some  constraints  on  models 
of  the  formation  and  evolution  of  the  outer  solar  system. 

6.1  The  MMT/ Megacam  Survey  for  small  KBOs 

I  have  devised  a  new  observational  method  which  allows  fast  photome¬ 
try  with  large  telescopes  with  standard  CCD  cameras.  I  am  able  to  achieve  high 
photometric  rates  (200  Hz)  on  tens  of  targets  simultaneously.  The  data  reduction 
techniques  for  this  kind  of  data  are  still  under  development.  The  amplitude  of  the 
noise  is  typically  larger  than  the  Poisson  noise,  and  it  displays  obvious  deviations 
from  normality.  However,  I  proved  this  method  is  suitable  for  gathering  a  large 
amount  of  precision  fast  photometric  data  in  few  observing  hours.  I  presented  a 
result  that  lowers  the  upper  limit  set  by  similar  occultation  surveys  by  more  than 
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one  order  of  magnitude  for  KBOs  D  >  1  km,  and  I  can  push  the  upper  limit  to 
D  >  700  m.  I  conhrmed  the  result  obtained  by  dedicated  Kuiper  belt  occultation 
surveys. 

The  high  speed  sampling  achieved  with  continuous  readout  mode  will  en¬ 
able  the  resolution  of  the  diffraction  features  of  any  candidate  events,  which  is  not 
possible  with  the  TAOS  project  due  to  the  lower  sampling  rate.  This  will  allow  us  to 
set  tight  constraints  on  the  physical  characteristics  of  an  occulting  system,  possibly 
breaking  the  degeneracy  between  impact  parameter,  size  and  distance  for  sub-km 
KBOs.  Furthermore,  continuous  readout  mode  enables  the  simultaneous  monitoring 
of  as  many  as  100  stars,  which  is  a  distinct  advantage  over  the  surveys  of  Roques 
et  al.  (2006)  and  Bickerton  et  ah  (2008),  where  only  two  stars  can  be  sampled  at  a 
time. 

This  observational  technique  has  proved  useful  in  testing  telescope  per¬ 
formance  and  addressing  guiding  issues  and  it  was  used  at  the  MMT  to  test  the 
drive  servos.  Furthermore  this  observational  method  is  a  promising  technique  for 
ground-based  high  precision  photometry  of  bright  sources  with  large  telescopes  as 
it  addresses  many  issues  typically  encountered  in  observing  bright  targets  (Gillon 
et  ah,  2008).  Saturation  is  avoided  without  resorting  to  refocusing,  it  involves  no 
overhead  due  to  readout  and  with  a  camera  like  Megacam,  with  a  large  held  of  view, 
it  allows  the  observation  of  many  stars  at  a  time,  guaranteeing  the  presence  of  a 
good  number  of  comparison  stars  that  can  be  used  to  achieve  high  precision  relative 
photometry. 

Further  improvements  in  SNR  might  be  achieved:  I  am  exploring  a  htting 
photometry  package  that  uses  the  Expectation-Minimization  algorithm,  treating 
each  row  as  a  mixture  of  Gaussians,  to  better  separate  the  contribution  from  different 
sources.  A  possible  way  to  address  the  contamination  due  to  unresolved  sources  is 
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to  subtract  the  contribution  from  known  unresolved  sources  (identified  from  the 
stare-mode  image,  see  Section  2.3.1)  using  the  trends  identihed  in  the  de-trending 
phase  (Section  2.3.1).  Another  possibility  is  to  de-trend  the  lightcurves  recursively, 
while  allowing  a  variable  phase  offset.  Finally,  it  shall  be  noticed  that  Megacam 
will  become  available  for  observations  at  the  Magellan  Clay  Telescope,  from  where 
our  target  helds,  at  the  intersection  of  the  galactic  and  ecliptic  plane,  could  be 
observed  at  a  higher  elevation.  This  would  help  reduce  the  noise  introduced  by  cross 
contamination  and  differential  image  motion,  as  the  atmospheric  effects  I  encounter 
observing  at  high  air  masses  would  be  reduced. 

6.2  The  TAOS  KBO  survey 

I  presented  an  analysis  of  3.75  years  of  TAOS  data,  comprising  5  x  10^  star- 
hours  of  three-telescope  4  or  5  Hz  time  series.  I  searched  for  occultations  of  our  star 
targets  by  KBOs  in  order  to  constrain  the  size  distribution  of  KBOs,  particularly  in 
the  1-to-lO  km  region,  which  is  currently  out  of  reach  of  direct  observation  surveys. 
I  described  the  TAOS  data  and  analyzed  the  dependence  of  the  efficiency  of  event 
recovery  on  various  parameters.  Occultations  near  opposition  lead  to  a  higher  event 
rate  for  TAOS,  even  after  taking  into  account  the  increased  recovery  efficiency  for 
small  objects  where  the  angle  from  opposition  is  larger  and  the  relative  velocity  of  the 
KBOs  is  lower.  The  largest  number  of  recoveries  is  obtained  for  a  stellar  magnitude 
TTtaos  ~  12.5.  Although  this  magnitude  is  close  to  the  limiting  magnitude  of  the 
survey,  the  abundance  of  Mtaos  ~  12.5  stars  in  the  TAOS  helds  compensates  for 
the  increased  SNR,  and  thus  detectability,  for  brighter  targets. 

I  found  no  occultation  events  in  our  data.  This  allowed  me  to  set  upper 
limits  to  the  number  density  of  KBOs  that  are  stringent  enough  to  be  compared 
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with  the  most  popular  models  for  the  formation  and  evolution  of  the  Kuiper  belt.  I 
considered  three  theoretical  models,  Pan  &  Sari  (2005,  hereinafter  PS05),  Kenyon  & 
Bromley  (2004,  hereinafter  KB04),  and  recent  simulations  by  Benavidez  &  Campo 
Bagatin  (2009,  hereinafter  BCB09),  all  of  which  describe  the  present  size  distribution 
of  the  Kuiper  belt  and  I  set  constraints  on  these  models.  This  is  the  hrst  comparison 
of  occultation  data  with  model  results. 

Our  result,  particularly  when  compared  with  PS05  and  KB04,  suggests  that 
the  Kuiper  belt  is  populated  by  fragile  bodies,  and  that  the  effect  of  the  migration 
of  Neptune  played  an  important  role  in  its  formation.  The  models  described  by 
BCB09  incorporate  specihc  fragmentation  prescriptions  (Farinella  et  ah,  1982;  Benz 
&  Asphaug,  1999)  but  do  not  include  a  model  for  Neptune’s  stirring.  The  BCB09 
models  depend  on  assumptions  on  the  initial  conditions:  the  location  of  the  hrst 
break  in  the  power  law  distribution  for  large  objects  and  the  slope  of  the  power  law. 
With  the  initial  conditions  used  by  the  authors  (Hb  =  100  km  and  gL  =  4.0)  TAOS 
allows  all  four  models  presented  in  BCB09. 

Using  a  simple  parametric  model,  where  the  size  distribution  is  described 
by  three  consecutive  power  laws,  and  hxing  the  slope  on  the  large  end  size  and  the  lo¬ 
cation  of  the  hrst  break  to  the  best  ht  from  direct  observations  (gL  =  4.25  km,  Ubi  = 
75  km),  and  the  slope  of  the  small  side  to  that  expected  from  a  collisionally  evolved 
strength-dominated  population  (gs  =  3.5)  I  can  constrain  together  the  intermediate 
slope  gi  and  the  location  of  the  second  break  /lb2,  as  shown  in  Figure  4.8.  As  direct 
surveys  are  currently  not  sensitive  to  KBOs  smaller  than  D  ~  28  km  the  TAOS 
occultation  data  can  for  the  very  hrst  time  probe  this  region  of  the  size  spectrum. 

I  also  considered  the  Jupiter  Family  Comets.  Assuming  the  JFCs  are 
injected  into  their  present  orbits  from  one  of  the  Kuiper  belt  populations.  Classical 
Belt,  Plutinos,  or  from  the  scattered  disk,  I  compared  the  upper  limit  derived  from 


109 


our  survey  to  the  estimates  of  the  number  of  objects  derived  using  the  number 
of  JFCs  respectively  by  Levison  &  Duncan  (1997)  for  a  population  of  cometary 
precursors  entirely  in  the  CB,  that  of  Morbidelli  (1997)  for  Plutinos  and  of  Volk  & 
Malhotra  (2008)  for  a  progenitor  population  in  the  scattered  disk.  I  can  rule  out  the 
CB  as  unique  precursor  family  if  the  progenitors  are  D  >  3  km.  This  preliminary 
result  conhrms  that  occultation  surveys  can  help  understanding  the  origin  of  JFCs. 

TAOS  has  operated  for  now  over  four  years  observing  continuously  with 
three,  and  now  four,  50  cm  aperture  telescopes.  TAOS  is  only  marginally  sensitive 
to  sub-km  KBO  occultations,  but  we  were  able  to  prove  that  the  low  sensitivity  at 
sub-km  sizes  is  more  than  compensated  for  by  the  vast  exposure  of  which  TAOS  is 
capable. 

6.3  The  TAOS  Sedna  and  Scattered  Disk  Survey 

The  TAOS  data  is  a  great  resource  for  probing  regions  of  the  outer  so¬ 
lar  system  that  are  too  distant  to  allow  direct  detections.  At  the  distance  of  the 
scattered  disk,  extended  scattered  disk,  and  Oort  clond,  objects  as  large  as  even 
thousands  of  kilometers  would  elude  direct  detection,  while  occnltation  have  the 
potential  of  reaching  these  remote  regions.  We  are  scanning  the  dataset  described 
in  Chapter  3  searching  for  shallow  flux  drops  affecting  several  consecntive  points, 
compatible  with  occultations  by  objects  in  orbit  between  100  and  1000  AU.  We 
detected  no  such  events  in  our  data. 

This  work  reproduces  what  was  presented  by  the  TAOS  collaboration  in 
Wang  et  ah  (2009)  with  two  major  improvements.  The  dataset  we  analyze  here 
is  larger  by  a  factor  of  ~  3  than  the  dataset  analyzed  in  Wang  et  ah  (2009).  We 
upgraded  the  search  algorithm  used  in  Wang  et  al.  (2009).  The  algorithm  we  use 
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here  accounts  for  non  stable  noise  characteristics  of  our  time  series  and  thus  it  allows 


us  to  exactly  calculate  the  probability  of  a  particular  conhguration  of  flux  values  in 
our  lightcurves.  Furthermore  we  expect  an  increase  in  recovery  efficiency.  At  this 
time  we  are  able  to  estimate  a  lower  limit  to  our  improvement  of  the  result  published 
in  Wang  et  ah  (2009). 

We  want  to  include  realistic  models  for  the  distribution  of  scattered  disk 
objects,  and  objects  in  Sedna-like  orbits  thus  possibly  setting  constraints  on  the 
genesis  of  the  orbit  of  Sedna  and  on  the  early  evolution  of  the  solar  system. 


Ill 
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